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ABSTRACT 


High  energy  laser  weapons  have  been  evolving  progressively  in  recent  years.  These 
weapons  deliver  high-intensity  beams  to  a  target  and  can  instantly  destroy  or  bum  it. 
They  may  cause  potential  threats  to  Navy  ships,  computer  networks,  guided  missiles,  and 
satellites  in  orbit.  In  order  to  reduce  our  military’s  vulnerability  to  high  energy  laser 
weapons,  one  possible  countermeasure  is  to  rotate  or  rock  the  object  itself  when  it  is  hit 
by  the  laser  beam. 

The  main  purpose  of  this  thesis  is  to  investigate  the  relationship  between  the 
speed  of  a  rotating/dithering  laser  beam  and  the  maximum  temperature  rise  induced  by 
the  laser  beam  on  a  finite  solid.  We  have  investigated  extensively  the  numerical 
solutions  for  the  transient  temperature  rise  in  both  one-dimensional  (1-D)  and  two- 
dimensional  (2-D)  finite  solids  due  to  rotating/dithering  laser  beams.  Our  mathematical 
approaches  include  the  eigenfunction  expansion  method,  the  Crank-Nicolson  method,  the 
Fast  Fourier  Transform  method,  and  COMSOL  for  1-D  and  2-D  cases.  We  have 
employed  COMSOL  to  solve  the  3-D  nonhomogeneous  heat  equation. 

This  thesis  provides  the  first  study  that  we  know  of  on  the  effect  of 
rotating/dithering  laser  beams  on  a  finite  target.  Our  results  are  consistent  with  previous 
analytical  studies  on  semi-infinite  regions.  The  quantitative  relationship  between 
maximum  temperature  rise  and  laser  beam  rotating  speed,  which  is  presented  in  this 
thesis,  can  be  used  as  a  general  guide  for  adjusting  the  speed  of  rotation  of  the  target  in 
order  to  prevent  the  maximum  temperature  rise  from  reaching  the  melting  point  of  the 
target. 
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I. 


INTRODUCTION 


Laser-induced  heating  has  been  widely  studied  [1,  2,  3,  4],  Most  of  these  studies 
are  restricted  to  a  moving/scanning  laser  beam.  Recently,  the  effect  of  a 
rotating/dithering  laser  beam  on  the  temperature  rise  of  a  semi-infinite  domain  has  been 
studied  [5],  However,  in  many  realistic  problems  the  geometry  of  an  object  is  finite.  So, 
in  this  thesis,  we  extend  the  previous  studies  to  a  finite  domain  and  pinpoint  the  effect  of 
a  rotating/dithering  laser  beam  on  the  temperature  rise  of  a  finite  object. 

The  laser  beam  can  be  considered  as  a  moving  disc  heat  source  with  a  Gaussian 
distribution  of  heat  intensity.  The  analytical  solutions  can  be  used  to  determine  the 
temperature  rise  distribution  in  and  around  the  laser  beam  heat  source  on  the  work 
surface  as  well  as  with  respect  to  depth  at  all  domains  close  to  the  heat  source.  The 
analytical  and  numerical  solutions  from  mathematical  approaches  provide  a  better 
appreciation  of  the  physical  relationships  between  the  relevant  laser  parameters  and  the 
thermal  properties  of  the  work  piece. 

The  main  purpose  of  this  thesis  is  to  investigate  the  relationship  between  the 
rotating/dithering  laser  beam  speed  and  the  maximum  temperature  rise  induced  by  the 
laser  beam  on  a  finite  solid.  In  Chapters  II  and  III,  we  present  the  analytical  and 
numerical  methods  for  obtaining  a  transient  temperature  distribution  with  various 
methods,  which  include  the  eigenfunction  expansion  method,  the  Crank-Nicolson  method 
and  the  Fast  Fourier  Transform  method  (FFT).  In  Chapters  III  and  IV,  numerical 
solutions  for  one-dimensional  (1-D)  and  two-dimensional  (2-D)  nonhomogeneous  heat 
equations  from  MATALB  and  COMSOL  are  given.  We  compare  different  methods  and 
investigate  the  relative  errors  of  numerical  solutions  obtained  in  MATLAB  and 
COMSOL.  In  Chapter  V,  we  have  employed  COMSOL  to  solve  the  three-dimensional 
(3-D)  problem. 

The  objective  of  this  thesis  is  to  provide  insights  on  developing  objects  of 
adequate  rotating  speed  against  direct  energy  weapons  (DEWs).  The  simulation  results, 
especially  the  quantitative  relationship  between  maximum  temperature  rise  and 


1 


rotating/dithering  speed,  can  be  used  as  a  blueprint  to  adjust  the  speed  of  rotation  of  the 
target  against  DEWs  to  prevent  the  temperature  rise  from  reaching  the  melting  point  of 
the  target. 
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II.  ANALYTICAL  SOLUTION  FOR  A  TRANSIENT 
TEMPERATURE  DISTRIBUTION  DUE  TO  A  ROTATING  OR 
DITHERING  LASER  BEAM  IN  A  FINITE  SOLID:  EIGEN¬ 
FUNCTION  EXPANSION  METHOD 


Consider  a  laser  beam  source  used  to  heat  a  1-D  finite  rod.  A  1-D  heat  equation 
can  be  written  as  [1] 


du  d2u  aT  ,  . 

—  =  aT — —  H — —q(x,t)  0  <x<L 

dt  T  dx2  K* 


(Eq.  2.1) 


where  u(x,t )  is  the  temperature  rise  with  respect  to  the  ambient  temperature,  q(x,t)  is 


the  heat  source  created  by  laser  beam,  aT  is  the  thermal  diffusivity  with  units 


m 


and 


Kt  is  the  thermal  conductivity  with  units 


W 
m  •  k 


Two  assumptions  are  made  in  this  analytical  approach.  First,  heat  losses  by 
radiation  are  negligible  as  compared  to  the  intensity  of  the  incident  laser  beam.  Second, 
thermal  properties,  such  as  thermal  conductivity  KT  and  thermal  diffusivity  aT  ,  are 

considered  constant  and  evaluated  at  an  average  temperature.  The  second  assumption 
makes  this  heat  equation  a  linear  problem. 

For  this  1-D  rod,  we  impose  insulating  boundary  conditions  (BCs) 

du  du 

— (0 ,t)  = — (Lx,t)  =  0  ,  which  assumes  that  no  energy  escapes  into  the  air  at  the 

dx  dx 

air/material  interface.  This  is  a  good  approximation  for  most  materials  under 
consideration,  as  heat  flow  by  conduction  through  the  material  far  exceeds  the  heat  loss 
by  radiation  or  convection  at  the  air/material  interface. 

The  initial  condition  (IC)  is  u(x,  0)  =  0 ,  which  reflects  the  fact  that  there  is  no 
temperature  rise  with  respect  to  the  ambient  temperature  before  the  laser  hits  the  rod. 


3 


The  heat  source  q(x,t)  is  considered  as  a  laser  beam  of  temporal  continuous  wave 
(CW)  and  spatially  modeled  by  a  Gaussian  distribution,  which  corresponds  to  the 
theoretical  TEM^  mode  of  the  laser.  The  term  TEM  comes  from  the  acronym 

“transverse  electromagnetic  mode.”  The  subscript  of  TEM  specifies  the  number  of  nodes 
generated  by  a  slight  misalignment  of  the  mirrors  located  in  the  laser  cavity  [1,6].  Figure 

1  shows  some  transverse  modes  and  the  simplest  modes,  TEM^  is  used  in  this  thesis; 
the  heat  source  generated  by  the  laser  beam  is  a  Gaussian  distribution. 
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Figure  1.  Different  transversal  modes  in  a  laser  spot  [From  1] 


For  a  dithering  laser  beam  in  TEMm  mode,  the  heat  source  q(x,t )  can  be  written 


as 


q(x,t)  =  I0e 


-- k-j(x-xc(t))2 


,  .  .  2nt  L 

,xc(t)  =  x0  +  a sin  —  ,  *0=y 
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Here,  xc(t)  is  the  position  of  the  dithering  Gaussian  beam  and  x0  is  the  initial 
position  of  the  laser  beam.  In  the  above  formula,  the  center  point  of  the  rod  is  picked.  r0 
is  the  effective  radius  of  the  laser  beam.  I0  is  the  intensity  of  the  laser  beam  at  the  center 

of  the  heat  spot  after  it  is  absorbed  by  the  material.  A:  is  a  constant  used  for  the  Gaussian 
model  [2],  Figure  2  depicts  a  dithering  laser  beam  shining  on  a  1-D  rod. 
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Figure  2.  Dithering  laser  beam  on  a  1-D  rod 


We  will  apply  the  eignefunction  expansion  method  to  solve  this  1-D  heat  equation. 
Generally  speaking,  the  method  of  eigenfunction  expansions  consists  of  building  up  the 
solution  of  a  boundary  value  problem  as  a  sum  of  eigenfunctions  of  a  Helmholtz  problem. 
To  begin  with,  consider  a  simple  1-D  mathematical  model  in  the  Cartesian  coordinate  for 
solving  the  linear  and  homogeneous  heat  equation  with  initial  conditions  (ICs)  and 
boundary  Conditions  (BCs)  described  above  where  no  heat  source  q(x,t)  is  involved,  the 
1-D  heat  equation  becomes: 


du 

dt 


=  u . 
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To  make  the  derivation  process  simpler,  set  otj  =  1 .  The  method  of  separation  of 

variables  (SOV)  is  applied  in  solving  this  homogeneous  heat  equation  [1],  We  will  move 
on  to  a  nonhomogeneous  solution  with  heat  source  q(x,t)  by  applying  the  method  of 
eigenfunction  expansion. 

It  is  intuitional  that  the  temperature  u  is  a  function  of  position  and  time, 
therefore  u  =  li(x,t) .  In  the  method  of  SOY,  we  attempt  to  determine  solutions  in  the 
product  form, 

u(x,t)  =  X(x)T(t)  (Eq.  2.2) 

where  X(x)  is  a  function  of  space  x  alone  and  T(/)  is  a  function  of  time  t  alone.  This 

SOV  reduces  a  partial  differential  equation  (PDE)  to  several  ordinary  differential 
equations  (ODEs). 

For  a  heat  equation  without  heat  source: 

^  =  V2 „(*,()  (Eq.  2.3) 

We  substitute  the  assumed  product  form  of  Eq.  2.2  into  this  heat  equation  and  get 
X(x)T\t)  =  X"{x)T(t) 

We  can  actually  “separate  variables”  by  dividing  both  sides  of  this  equation  by 
X(x)T(t)  to  obtain 


r  _  x” 

T  ~  X 


(Eq.  2.4) 


Now  the  variables  have  been  “separated”  in  the  sense  that  the  left-hand  side  of  the 
equation  is  only  a  function  of  timet  and  the  right-hand  side  function  is  a  function  of 
space  x .  Since  the  variables  t  and  x  are  independent  of  each  other,  the  only  way  to  get 
equality  is  to  have  the  functions  on  both  sides  of  (Eq.  2.4)  constant  and  equal.  Thus, 


T' 

T 


V" 

X 


=  -A 
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where  A  is  an  arbitrary  constant  called  the  separation  constant.  Later  we  will  explain 
why  this  separation  constant  is  a  negative  value  with  A>0 .  Since  we  treat  the  laser  beam 
as  a  heat  source,  and  there  is  no  heat  source  involved  yet  in  solving  the  homogeneous 
equation,  it  is  intuitive  that  the  temperature  decreases  as  a  function  of  time,  so 

T'(t) 

— ^  =  -A<0  (Eq.  2.5) 

T(t ) 

(Eq.  2.5)  is  a  first-order  linear  homogeneous  differential  equation  with  constant 

yf 

coefficient.  This  ODE  can  be  solved  directly  by  seeking  exponential  solutions,  T  =  e  . 
By  substitution,  the  characteristic  polynomial  is  r  = -A.  Therefore,  the  general  solution 
of  (Eq.  2.5)  is 

T(t)  =  ce~ ^  (Eq.  2.6) 

Where  c  is  an  arbitrary  constant.  The  time-dependent  solution  is  a  simple 
exponential.  Recall  that  A  is  the  separation  constant,  which  for  the  moment  is  arbitrary; 
however,  eventually  we  will  find  out  that  only  certain  values  of  A  are  allowable. 
If  A  >  0 ,  the  solution  exponentially  decays  as  t  increases  because  of  the  negative  sign; 
if  A  <0 ,  the  solution  exponentially  increases;  if  A  =  0 ,  the  solution  remains  constant  in 
time.  Therefore,  A  >  0  gives  a  reasonable  solution  since  temperature  decreases  as  time 
goes  by. 

Our  next  move  is  to  consider  the  right-hand  side  of  (Eq.  2.4): 

X—  =  -A  (A  >  0)  (Eq.  2.7) 

This  is  a  second-order,  constant  coefficient  homogeneous  ODE,  substituting 
X  ( x )  =  erx  into  (Eq.  2.7)  yields  the  characteristic  equation  r  +  A  =  0 . 
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Since  X  >  0 ,  exponential  solutions  have  imaginary  exponents,  X(x)  =  e±l^x .  In 

iO 

this  case,  the  solution  oscillates.  Recalling  Euler’s  formula,  elu  —  cos(0)  +  i  sin(i9) ,  so 
the  choices  COS^VXxj  and  sin(VXx)  are  made  to  get  the  desired  real  number  solutions. 
The  general  solution  for  (Eq.  2.7)  is: 

X(x)  =  Acos(\[Jx)  +  Bsm(JJx)  (Eq.  2.8) 

This  is  an  arbitrary  linear  combination  of  two  independent  solutions  where  A  and 
B  are  just  arbitrary  constants. 

Both  ends  of  the  1-D  rod  are  insulated;  the  BCs  on  x  direction  states: 

du  A  ,  du  A 

—  -0  and  —  =0 

dX  x=0  dx  x=Lx 

This  is  equivalent  to 

X'(0)  =  0  and  X\Lx)  =  0  (BCs) 

The  derivative  of  (Eq.  2.8)  is 

X  '(x)  =  -AU  sin(VXx)  +  B\[X  cos(VXx)  (Eq.  2.9) 

We  can  plug  the  first  BC  X  '(0)  =  0  into  (Eq.  2.9)  to  obtain  the  solution 
X  '(0)  =  0  which  implies  that  B=0.  SoX(x)  is  a  function  of  cosine  only  since  the  sine 

term  vanishes.  Another  BC  implies  that  A  sin(VXZx )  =  0  .  To  avoid  the  trivial 

solution  X (0)  =  0  ,  we  take  A  =  1  ,  which  forces  sin(y/ALx )  =  0  .  Since  the  sin 
function  vanishes  at  the  integer  multiples  of  n ,  we  conclude  that 

=  n  —  0,±1,±2, . 

and  so 
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rYlTlX, 


X  =  Xn  (x)  =  cos(-^— ),  n  =  0,1,2,. 


Note  that  for  negative  values  of  n ,  we  obtain  the  same  solution;  hence,  solutions 
corresponding  to  negative  n's may  be  discarded  without  loss.  These  Xn (x)  are  from  the 
eigenfunction  of  the  1-D  heat  equation.  We  now  go  back  to  (Eq.  2.6)  and  substitute 


X  = 


f  \2 

nn 

\Lx  j 


to  get 


H7T 


T„(t)  =  cne 


v4v 


Now  we  return  to  (Eq.  2.1).  We  solve  this  problem  using  the  method  of 
eigenfunction  expansion,  and  hence  start  by  assuming  that  u  has  an  expansion  in  terms  of 
the  eigenfunctions  Xn  (x) .  Thus, 


,  s  x- '  /  x  .TITCX . 

O,0  =  2^(0c  os(— — ) 


u 


(Eq.  2.10) 


eigenfunction  exp  ansion 

where  the  coefficients  un  (t)  are  functions  of  t .  We  now  expand  q  as 

q(x,t)  =  YJq,M)xn(x) 


(Eq.  2.1 1) 


After  substituting  (Eq.  2.10)  and  (Eq.  2.11)  into  the  1-D  inhomogeneous  heat 
equation  (Eq.  2.1),  we  obtain 


2 ui  4)4 (x)  =  aTY,u„(t) 


f  \ 

nn 

2 

v  ) 

K 


T  n 


This  yields  the  differential  equation  of  the  coefficients  un  (t) , 


un  \t)  =  aT 


f  \ 

2 

nn 

L 

\  X  J 

CCrp 


«,(0+Tf^(0 

A  , 


(Eq.  2.12) 
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The  initial  condition  for  this  equation  is  un  (t  =  0)  =  0 . 


Now  the  problem  is  to  find  qn(t) .  We  can  use  the  property  of  orthogonality  of 


M7TX, 


eigenfunctions.  By  multiplying  eigenfunction  Xn  (x)  which  is  COS(— — )  and  then 

x 

integrating  on  both  sides  of  Eq.  2.1 1,  we  get 


|  q(x,  t)Xn  (x)dx  =  \YJcln  (t)xn  ( X )  •  Xn  (x)dx  =  qn  (t)  |  Xn 2  (x)dx 


0  « 


,  .2  nnx. 

£  l  +  c°s(— — ) 


L-  dx  =  qn{t)^z 


2  Lr  -4(— c(0)2 

1  -  r0 


,nnx^ 


Therefore,  qn  (0  =  —  J  q(x,  t)Xn  (x)dx  =  —  jl0e  r°~  cos (^j^dx 


X  0 


X  0 


To  obtain  qn{t )  is  very  expansive  because  the  x-xc(t)  term  makes  the  integration 
time-dependent.  Here  we  use  change  of  variables  to  make  this  calculation  easier, 

X  =  x-xc(t) 
dx  =  dx 

thus, 


-  Lx-Xc(t)  k  -2 

2  f  r2x  ,tl7T 


In  (0  =  —  J  he  r° 2  cos(—  (x  +  xc  (t)))d X 


x  -Xc(t) 


We  introduce  a  constant  d ,  which  is  the  effective  radius  of  the  laser  beam,  and 

k  ~2 
~—X 

d  is  about  3  to  5  times  of  the  standard  deviation  S  of  I0e  0  .  Instead  of  taking  integral 

from  x  =  -xc(t)  to  Lx-xc(t)  ,  using  the  assumption  above  to  change  the  integral 

from  x  =  -d  to  x  =  d ,  since  the  laser  beam  is  a  very  concentrated  Gaussian  distribution, 
so  taking  integral  with  respect  to  a  large  range  is  unnecessary;  we  only  need  to  integrate 
within  the  effective  laser  range  which  is  from  -d  to  d .  Figure  3  shows  the  effective 
radius  of  the  laser  beam. 
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I, 


-d 
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d 


Figure  3.  Effective  radius  of  laser  beam,  d  as  a  multiple  of  standard  deviation  8 

We  also  recall  that  the  trigonometry  formula  of  cosine  in  sum  and  difference 
formula, 

cos(w  ±  v)  =  cos  (u)  cos(v)  +  sin(w)  sin(v) 
then  qn  (t)  becomes 


The  integral  in  the  second  term  is  zero  since  the  integrand  function  is  an  odd 
function  and  the  integration  interval  is[-d,  d].  Thus,  qn(t)  simplifies  to 


Once  qn  (/)  is  found,  we  have  to  compute  un  (t)  using  the  method  of  integrating 

factor. 

Recall  the  method  of  integrating  factor  in  solving  first-order  ODE: 
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yXt)  +  a-y(t)  =  r(t) 

Here  a  is  a  constant  and  r(t)  is  a  function  of  time.  The  first  step  of  integrating 


at 

factor  is  to  multiply  e  on  both  sides  of  this  equation: 
eat  -[yW  +  a-yd^e01  -r(t) 


This  equation  is  equivalent  to 


We  can  integrate  on  both  sides  of  the  equation  and  get 

t 

eat  •  y{t)  =  7(0)  +  J  eas  •  r(s)  ds 

t- 0 

where  s  is  a  dummy  variable.  We  multiply  e  to  obtain 

y(t)  =  .y(O)  •  e~at  +  e~at  •  J  eat  •  r(t)  dt 

t= 0 

Now,  we  can  compute  un  (t)  of  Eq.  2.12  use  the  integrating  factor  method  to  get 
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where  the  initial  condition  w„(0)  has  been  applied. 


Once  un  (t)  is  calculated,  we  can 


find  u(x,t)  =  =  2«.(')cos(— ) 

«  «  Jt 


After  we  conduct  several  experiments  in  MATLAB,  we  find  that  the  direct 
summation  is  very  slow;  the  computational  cost  of  this  eigenfunction  expansion  method 
is  extremely  high  in  a  3-D  nonhomogeneous  heat  equation  problem.  An  alternative  way 
to  compute  the  series  is  to  take  advantage  of  the  Fast  Fourier  Transform  (FFT).  In  order 
to  do  so,  we  first  change  the  cosine  form  to  a  combination  of  exponential  forms. 


.rmx .  e 
cos{——)  =  — 


.  nnx  .  nnx 

11T  ,  1~l~ 


+  e 


Then  we  even  extend  u(x,t)  from  the  domain  (0 ,LX)  to  (~Lx,0)  and  make  it 
periodic  with  period  2 Lx  .  Figure  4  shows  how  to  expand  the  domain  in  FFT  method. 


-A  0  l. 


Figure  4.  Domain  expansion:  cosine  expansion  to  FFT 
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So,  we  rewrite  u  as 


u(x,t)  =  YJU„(t) 


f  .nnx 
1  Lv 


n= 0 


+  e 


.nnx  A 

1  u 


®  i - 

=  2_j  c„  (0 ' e  4  where 


c„(0  =  i 


K(0 
2  ’ 
H-Jt) 


n  >  0 
n<0 


and  then  we  apply  FFT.  We  will  provide  more  details  on  how  to  compute  the  solution 
with  FFT  in  Chapter  III. 
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III.  NUMERICAL  SOLUTION  FOR  A  TRANSIENT,  ONE¬ 
DIMENSIONAL  TEMPERATURE  DISTRIBUTION  IN  A  FINITE 
ROD  DUE  TO  A  DITHERING  LASER  BEAM 


In  this  section,  we  review  some  numerical  methods  for  solving  the  1-D 
nonhomogeneous  heat  equation. 

A.  THE  CRANK-NICOLSON  METHOD 


The  Crank-Nicolson  method  is  an  implicit,  unconditionally  stable  numerical 
method  used  to  solve  the  heat  equation.  It  is  second-order  accurate  in  time  and  second- 
order  accurate  in  space.  In  1947,  Crank  and  Nicolson  suggested  an  alternative  way  to 
utilize  centered  differences  scheme  [7].  To  illustrate  the  Crank-Nicolson  scheme,  we 
start  with  centered  difference  and  return  to  the  Taylor  series  expansion  for  / (x0  +  Ax)  and 

f(x0  -  Ax) . 


/(,+At)=/w+|/w+(Af/.w+M/P)W+... 

/(*  -  a*) = m  -  77  /  '<  v + ^-  / n  w  -  ^-/®  w + ... 

Subtracting  (Eq.  3.1)  from  (Eq.  3.2)  yields 

fix  +  Ax)  -  f(x  -  Ax)  =  2  Ax  /  '(x)  + (  Ax)3  /(3)  (^)  +  •  •  • 


We  can  get  f\x)  by  reorganizing  this  equation 


2Ax  6 


=0(  Ax)1 


This  leads  to  the  first-order  centered  difference  approximation 
f(x  +  Ax)  —  f{x  —  Ax) 


f\x) 


2Ax 


(Eq.  3.1) 

(Eq.  3.2) 
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This  approximation  is  with  truncation  error  0(  Ax)2,  an  order  of  (Ax)2  and  the 
approximation  is  accurate  and  consistent;  the  truncation  error  vanishes  as  Ax  — >  0 . 

By  adding  Eq.  3.1  and  Eq.  3.2,  we  have 

/(x  +  Ax)  +  /(x  -  Ax)  =  2  fix)  +  (Ax)2  /  "(x)  +  ( Ax)4  /(4)  (x)  + ..  . 

Again,  we  can  get  the  second-order  derivative  f'\x)  by  reorganizing  this  equation 


/  "(x)  =  /(*  +  M  Z  +  /(*  Z  **)  _  y  (4 )  +  _ 


(Ax)2 
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O(Ax)2 


and  this  leads  to  the  second-order  centered  difference  approximation 
f(x  +  Ax)  -  2/(x)  +  fix  -  Ax) 


/"(*)' 


(Ax)2 


After  we  understand  the  second-order  difference  approximation,  we  can  move  on 
to  the  Crank-Nicolson  scheme.  The  forward  difference  in  time  of  temperature 
T  =  T(x,  t )  can  be  written  as: 

dT  _  Jjt  +  At)-Tjt) 
dt  At 

The  Crank-Nicolson  method  may  be  interpreted  as  the  centered  difference 

around  t  +  — .  The  error  in  approximating  —  (t  +  — )  is  0(A?)2  [7],  Thus,  we  discretize 
2  dt  2 

the  second  derivative  at  t  +— with  a  centered  difference  scheme.  Since  this  involves 

2 

functions  evaluated  at  this  in-between  time,  we  take  the  average  at  t  and  t  +  At.  This 

dT  d2T 

yields  the  Crank-Nicolson  scheme  for  1-D  homogeneous  heat  equation —  =  aT  — - , 

dt  fix 


rj-ik  +  \  rj-ik 

At 


(/. , 


rpk  _L  T7^  rpk+ 1  rjrpk+l  .  rpk+l 


Ax 


Ax2 


(Eq.  3.3) 
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We  introduce  the  shorthand  notation  T(xi,tk)  =  Tik ,  which  is  a  matrix  notation 

used  to  discretize  x  and  t ,  because  the  Crank-Nicolson  scheme  involves  six  stencil 
points  rather  than  a  simpler  stable  method,  and  three  of  which  are  at  the  advanced  time, 
as  shown  in  Figure  5.  We  cannot  directly  march  forward  in  time.  Instead,  we  are  left 

with  a  traditional  system  of  equations  of  the  form  Ax  =  b  ,  where  x  is  the  solution  at  the 
(&  +  l),st  time  level  and  b  depends  on  the  solution  at  the  k-th  time  level.  Of  course,  it  is 
not  just  any  A ,  the  matrix  we  have  has  a  very  special  structure  which  is  to  be  shown  later. 

t  +  A  t 


t 

x-Ax  x  x  +  Ax 

Figure  5.  Implicit  Crank-Nicolson  scheme  [From  7], 

Now  consider  a  1-D  nonhomogeneous  heat  equation  with  a  heat  source  / (x,t) : 

—  =  ^  +  f(x,t)  ,  0  <  x  <  1 
dt  dxlJ 

The  BCs  are  insulated  at  both  ends  of  the  rod: 

8T_  _dT_  _Q 

dx  **  ~  dx 

The  IC  assumes  that  there  is  no  temperature  rise  with  respect  to  the  ambient 
temperature: 

r(x,0)  =  0 
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In  this  thesis,  we  study  the  case  where  the  heat  source  is  a  dithering  1-D  Gaussian  beam: 

-(x-xc  (t))2 

f(x,t)  =  I0e  r°  ,xc(t)  =  x0+asm^,  x0=^- 

where  /0  is  the  intensity  of  the  laser  beam.  As  before,  we  discretize  this  1-D  rod  into  n 
parts,  as  depicted  in  Figure  6: 


Xi=(i-^)Ax, 


Ax 


(i-o) 


i  =  1,2, 


Cr . O - o- 

0  X1 


F-i=Ti 


Figure  6.  Discretization  of  the  1-D  rod 


-o - o . o 

xn  i  xn+1 


T„=T 


n+1 


Here  in  the  Crank-Nicolson  scheme  and  later  in  the  matrix  notation,  the  BCs  are 
extended  into  ghost  grid  points  x_x  and  xn+l .  However,  in  this  finite  difference  method, 

the  BCs  link  temperature  at  x_x  to  the  temperature  at  x,  and  the  temperature  at  xn+1  to  the 
temperature  at  xn .  More  specifically,  the  BCs  are 


dr_\  .  r(xi,c)-r(x_i,Q  _  k_ 

S  *=0  ~  .  J  —  / 

ox  Ax 


dT_ 
fix  1-1 


r(x„+1,4)-r(x„,c) 

Ax 


=  o=>r  =t 

yJ^1n+ 1  1n 


for  all  k>\ 


Applying  the  Crank-Nicolson  scheme  to  the  heat  equation  with  a  heat 
source  / (x,  t ) ,  we  obtain 


rj~iK+l  rj-ik 

At 


a 

~2 


rpk  .  rjik 

1i+\~zli  +  1 

Ax2 


rj-ik  + 1  ^rj-i/t+L  ^  rj- r, 


+ 


k+ 1 


Lz+1 


nk+l 

i- 1 


Ax2 


+  f[f(xntk)  +  f(xi’tk+i)\ 
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tk  =  kAt 


i  =  \,2,....n-\,n  k  =  0,1,2, . 

Here,  {7^  j  are  known,  |7’/l+l  j  are  to  be  found.  We  start  with  a  simple  case  i  =  1  and 
rearrange  this  equation: 


rnk-Y  1  rpk  .  &  At 

1\  ~  1\  “l  I” 


]-ik  rj-iK  +  i  r^rjik+X  ^  rjik  + 1 


Ax2 


Ax 


+  y  [/(VJ  +  ZCMm)] 


where  T\  =  T*  and  7^+1  =  Tk+]  because  of  the  BCs.  After  changing  some  orders,  we  put 
7^+1  on  the  left-hand  side  of  the  equation  and  Tk  on  the  right-hand  side;  we  get 


*+ 1  ate  _t+1  ,  aAt  _t+1  _t  _  aAt  t  aAt  At 


fyr/t+l  ^  ^  rpK+\  rp/i  +  l  rri  k  .  ^  ^  rji  k 

1\  _  ->  A  2  +  7  2  4  —  M  +77TJ2  _  ~  .  2 

2Ax  2Ax  2Ax  2Ax 


T\  +—[f(xi’tk)  +  f(xi’tk+ 1)] 


To  simplify,  we  get 


^  aAt  ^k+t  ^  aAt  (.  aAt  „t~\  aAt  At 


1  +  - 


2Ax 


2  ±  1 


2Ax 


Tk+l  = 


2  2 


1  — 


2Ax2 


+— [/(*„4)+/(*,>4.,)]  (£«■  i) 


2Ax 


To  put  this  in  matrix  form,  we  get 


1  + 


aAt  aAt 


2  Ax2  2  Ax2 


aAt  aAt 


Jlx2 


ik+ 1 


■7/t  +  l 


2x1 


2  Ax  2  Ax 


rjik 

4 


1x2  LJ2  J2xl 


+  -y  [/(xi  >  4 )  +  /(xi > 4+i )] 


Likewise,  when  t  =  2 


r2  -r2*  _a 
At  _  2 


Tk -2Tk +  Tk  Tk+l 
Ax2  + 


-  2Tk+l  +Tk+l 
~Ax2 


+  ^[/(X2>4)  +  /(X2>4+l)] 


After  changing  some  orders,  we  put  T  on  the  left-hand  side  of  the  equation  and 
Tk  on  the  right-hand  side 
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ato  t*+i  (.  a&Tk+ iA 
'  Ax2  2 


2AX2  1 

:  —  Tk  + 
2Ax?  1 


c^t_Tk+ 1 


v  ^  y 

^  /vA  #  ^ 


2Ac‘ 


.2  3 


1 


CtlSl 


Ac2  J  2  2  Ac2 


eiAt  ^  At 


(^•2) 


^2  +  —  [f(X2’tk)+f(X2’tk+l)] 


In  matrix  form, 


aA t  ,  aAt  «At 
1  +  - 


2  Ax2  Ax2  2  Ax2 


aAt  aAt  aAt 


1x3 


n£+l 


n£+l 


p&+1 


2  Ax2  Ax2  2  Ax2 


1x3 


3x1 


At, 


+  —  [f(X2>tk)  +  f(X2>tk+l)] 


J3xl 


When  i  =  n  —  1 


Tk+l  _Tk  a 

1  JB-1  _ 

At  2 


rrik  r^rrik  .  rrik  rrik+l  r\rrtk+ 1  .  rrrk+1 

1n  ~  A1n-\  +  _ZL_  ~Z1n- 1  ~h  1  n-2 


Ax 2 


Ax2 


+  ^  [/(x«-i > 4  )  +  /(x«-i » 4+i )] 


Once  again,  we  put  +1  on  the  left-hand  side  of  the  equation  and  Tk  on  the  right 
hand  side 


oAt  ym 

2Ax2  "-2 


*4  oAt^l  oAt 
Ar2  -1 


oAt  _* 

: - T  + 

2A*2  "-2 


1- 


rj~k+ 


y 


oAt 

A? 


2A*2  " 
aSt  At 


(£fcw-l) 


4f-!  +^2  ^  +7  l/^Vl  ’  4  )  +/(XH-1  ’  4+1 )] 


In  matrix  form 


aSt 

2Ax2 


aAt 

2Ax2 


1  + 


«At  «At 


Ax2 


2Ax 


1x3 


£+1 

■2 

'^+1 
n-\ 
<k+ 1 


« 


3x1 


1- 


aAt  aAt 


Ax2  2  Ax2 


-11x3 


rrik 
1n- 2 


T 

n 

T 


k 

n- 1 
k 


+f  [/(Vo4)+/(Vo4+i)] 


J3xl 
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Finally,  when  i  =  n 


rrik+l  rrik 

_n _ 1  n 

At 


a 


rrik  rsrpk  .  rrik  rpk+l  <jrpk+\  .  rpk+l 

1n+\~Z1n  ~Yin- 1  j  1  n+ 1  ~Z1n  ^  1  n- 1 


Ax2 


Ax2 


+  —  \f(Xn’tk)  +  f(Xn’tk+ 1)] 


Here,  Tk+l  =  Tk  and  =  Tk+]  because  of  the  BCs. 

Putting  r/,+l  on  the  left  hand  side  of  the  equation  and  Tk  on  the  right-hand  side 
leads  to 


Q&  ^+i 

~2Ar  -1 


1+ 


oAt 


rjJz+ 


2Av^  " 


aAt  ^ 

: - 7*  + 

2A*^ 


1- 


oAt 

2A? 


Tt  +^[f(A„tk)+f(xn,tk+])]  ( Eq.n ) 


In  matrix  form, 


aAt 

2Ax2 


1  + 


aAt 

2Ax2 


Jlx2 


^it+i 
/2-1 
t£+ 1 


/2 


2x1  L 


aAt 

2Ax2 


1- 


aAt 

2Ax2 


Jlx2 


/7-1 


/2 


2x1 


+  [  f(Xn  >  4  )  +  /(Xn  >  4+1 )] 


Putting  all  these  n  linear  equations  Eq.  1,  Eq.  2, ...  Eq.  n- 1,  Eq.  n  together  in  matrix 
form,  we  arrive  at 

— /t+l  —k  — 

AnxnTn*'  =BnxJn*  + /„xl  =  6«xl 


where 


1 

+ 

_ 1 

rpk 

1\ 

/Oi,4)  +  /Op4+i) 

— t+l 

rpk+l 

12 

—k 

tr  - 

—  At 

/(X2>4)  +  /(X2>4+l) 

T  = 

If  ••• 

T  = 

,1  ... 

/  =  — 
2 

/(x«-o4)  +  /(x«-i,4+i) 

1 

V-i 

+ 

1 _ 

nx  1 

_ i 

nxl 

.  /(x„>4)  +  /(X>4+i)  _ 
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mi 


0 


0 


a/St  a/St 
2A?  2A? 


a/St  a/St 
_A?  2A? 


aSt  y+c<^ 

2/S?  Ax2  2/S? 

_  a/St  ^  a/St 

0 - -  1+ — - 

2A*2  2A*2 


6i2V  |  oA?  a/St 
2A<?  Ax2  2A? 


0  0 


2A*2  2A*2 


,  ccSt  t  cxSl  k  St r  y , y  .  .-I 

2A?J^  +2A?  !  +T[/<X|’'*)+/(J;’,w)] 

+f1_7^V  +f  [/fe4)+/fe4.,)] 

y  AX  J  ZlsX  z 


+(l +f  [/(W*)+/(W*,)] 

exA/t  k  f .  exAt  ^  ^  r  /v  .  /,/  n.  "I 

2A?  +  2A?  "  +T[/(X”’(‘)+/(X”’'“)] 
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Note  that  A  and  B  are  triangular  matrices.  We  can  use  the  sparse  matrix  solver  in 
MATLAB  to  create  them  and  they  are  easily  to  be  computed.  In  MATLAB,  we  simply 

—k+ 1 

use  T  =  A  \  b  to  find  the  temperature  rise  with  respect  to  time. 

Let  us  recall  from  Figure  2,  the  idea  of  a  dithering  laser  beam  heat  source 
/ (x,t)  on  a  1-D  rod  of  Eq.  2.1  can  be  written  differently: 


dT 


=  aT 
dt  dx 


d2T  a, 


+  ^f(x,t),0<x<Lx 


K, 


f(x,t)  = 


yflftcl2 


-(x-xc(t)Y 
,  2d2 


,  xc  (t)  =  x0  +  a  ■  sin 


^  2  7Tt  ^ 
period. 


,  x. 


_  4 


(Eq.  3.4) 


Where  d  is  the  standard  deviation  of  the  Gaussian  heat  source.  Figure  7  demonstrates 
what  the  final  temperature  rise  looks  like  for  time  =  1  using  the  Crank-Nicolson  method 
with  number  of  points  n=  200.  All  the  input  parameters  are  listed  in  Table  1.  The 
MATLAB  code  is  attached  in  Appendix  A. 


Input 

Value 

Unit 

aT 

1 

mA2/s 

Kt 

1 

W/(m*k) 

h 

1 

W/mA3 

d 

0.02 

m 

4 

1 

m 

x0 

0.5 

m 

a 

0.25 

m 

period 

1 

s 

Table  1.  1-D  rod  MATLAB  input  parameters  for  dithering  laser  beam 
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1-D  Crank-Nicolson  method  with  n=200  @  t=1 


Figure  7.  1-D  Crank-Nicolson  method 


B.  THE  FAST  FOURIER  TRANSFORM  (FFT)  METHOD 

Recall  that  in  (Eq.  2.10)  we  have  shown  that  the  solution  of  the  1-D 
nonhomogeneous  heat  equation  with  insulating  boundary  conditions  can  be  expressed  as 
eigenfunction  expansion,  where  the  eigenfunctions  are  cosine  functions. 

Now  we  show  how  to  use  FFT  to  implement  the  eigenfunction  expansion  for  the 
purpose  of  fast  summation.  This  can  be  achieved  in  several  steps. 

First,  we  even  extend  the  initial  condition  from  [0,1]  to  [1,2]  and  then  make  it 
periodic  with  period  2.  After  that,  we  map  the  interval  [0,2]  to  [0,2;r]  .  By  doing  this,  the 
Fourier  series  will  only  contain  cosine  terms. 
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Second,  we  apply  the  discrete  Fourier  transform  to  the  heat 

.  du  d2u  ... 
equation  —  =  —  +  f(x,t): 
at  ox 


F 

du 

=  F 

~  d2u~ 

_  dt  _ 

dx2 

+  F[f(x,t)] 


where  F  denotes  the  Fourier  transform  operator. 

Now  we  need  to  calculate  Fourier  transformations  of  derivatives  of  u(x,t) .  We 
begin  by  recalling  the  spatial  Fourier  transform  of  u(x,t) : 

1  2?r 

F\u(x,t)\  =  u{k,t)  =  —  f  u(x,t)elk7UXdx 

?77-  J 


Note  that  this  is  also  a  function  of  time;  it  is  an  ordinary  Fourier  transform  with 
time  t  fixed.  To  obtain  a  Fourier  transform  in  space,  we  multiply  elkx  and  integrate. 
Spatial  Fourier  transforms  of  time  derivatives  can  be  derived  easily  because  the  spatial 
Fourier  transforms  of  a  time  derivative  equals  the  time  derivative  of  the  Fourier 
transform: 


F 

du(x,  t ) 

1  J  du(x,t)  cikxxdx  _  d 

i  2^ 

—  f  u(x,tyknxdx 

dt 

2k  l  dt  dt 

2tt  J 

L  o 

d_ 

dt 


u(k,t ) 


For  spatial  Fourier  transform  of  spatial  derivatives,  the  method  of  integrating  by 
parts  can  be  used: 


F 


du(x,  t ) 
dx 


1  J  du(x,t)  eik„xdx  =  ujxj)e^ 
2 k  *  dx  2k 


-ikK 


2tl 

—  J  u(x,t)e,kKX 
2k  n 


dx 


=  -ikKu(k,t) 


Here,  we  assume  w— »0asx— »0or2;r,  and  then  the  first  term  of  this  equation 
vanishes. 
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In  general,  the  Fourier  transform  of  the  nth  derivatives  of  a  function  with  respect 
to  x  equals  (- ik 7t)n  times  the  Fourier  transform  of  the  function,  assuming  that 
u(x,t)  — >  0  sufficiently  fast  as  x  approaches  periodic  endpoint  x  — >  0  and  x  — >  2n . 


Likewise,  Fourier  transforms  of  a  second  derivative  can  be  obtained: 

=  (-/Avr)2  u(k,t)  =  -k27r2u(k,t) 


F 

d2u(x,  t ) 

=  -iknF 

du(x,  t) 

d2x 

8x 

dw  A  A  A  A 

Finally,  we  can  conclude  that  —  =  {-ikn)1  u  +  f(x,t)  =  -k27r2u  +  f(k,t)  .  The 

dt 

factors  here  is  due  to  the  fact  that  we  map  the  interval  [0, 2]  to  [0,  2tt]  .  Using  the 
method  of  integrating  factor  and  multiplying  ek  K  1  on  this  equation,  we  find 

k2n2t  du  .  2  2  jdn2t*  .  jdn2t  2-/7  j.\ 

e  —  =  —k  7i  e  u  +  e  f(k,t) 
dt 

We  rewrite  this  ODE  as 


j^ueM)  =  eMf(k,t) 


Integrating  both  sides  yields 


u(k,t)  =  u(k,  0)<f*V'  +e-kWt\ek2^sf(k,s)ds  =u(k,0)e-k2*2t  +  J/V(s'-°/(A:,v)c/v 

o  o 

Here,  we  make  the  change  of  variables  t  =t-s  to  conduct  our  following  derivation. 


This  equation  becomes 

t 

i(k,t )  =  u(k,0)e~k  n  ‘  +  Je~k  n  T  f(k,t  —  T)dr 


m 


In  order  to  estimate  this  integral,  we  approximate  f{k,t-  r)  by  a  linear  function: 
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Then,  it  follows  that 


i 

u(k,  t )  =  u(k,  OK* V'  +  J  e'*Vr  f(k,  t-T)dr 


i 

\(k,  0)eM  +  j  e 


-k  n  x 


f(k,t)  + 


f(k,0)-f(k,t ) 


dr 


=  u(k,  OK  + 


Z-n.  +\  1  /i  ,  f(k,0)-f(k,t)  \ 

)+ — 7 — r 


rc/r 


=  u(*,  OK*V?  +  /(*,  0  ) 

K.  7T  '  ' 


+ 


f(k,0)- f(k,t) 


k V 

t  -k2^2!  1 

e  + 


Ar2;r2 


fcV 


«:  n  '  ’ 


+ 


f(k,0)-f(k,t) 


1  1  +  ffcV  _*v, 


A:V4  A:4;r4 


i-(mVKV( 


*V 


Once  w(A:,t)  is  found,  apply  the  inverse  discrete  Fourier  transform  to  it  to 
get  u(x,t) . 

In  our  MAT  ALB  code,  as  attached  in  Appendix  B,  t  is  replaced  by  dt  because  in 
each  step  along  the  time  direction  we  march  dt  .  The  first  component  in  or 

^corresponds  to  the  case  where  A:  =  0.  We  apply  H’Lospital’s  rule  to  find  the  value 
when  k  =  0 .  More  specifically,  we  have  used  the  formulas: 


lim- 

k^°  kZ7T 


1  /  9  7  \  A=k 


.¥,!"■  1-e 


-At 


r  0  ^ 


lim: 

i->0  T 


-type 


=  lim- 


<-0* 


-At 


1  — (1 j4Vff  l-(l  +  q)e 

™  A:4^-4  ™  T2 


y 


a-»o  ) 


=  t 


-type 


te  1  +  t(l  +  tX)e 


y 


=  lim- 

A->0 


-a* 


2/1 


y0  a 

-type 

y 
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Temperature 


=  X-  lim  \t2eXt  + 1 \telt  - 1(  1  +  tA)e~Xt  ]} 
=  lim  |  t2e~Xt  + 1  (-At2 


Figure  8  shows  the  result  using  FFT  method  in  MATLAB  where  all  the  input 
parameters  are  from  Table  1. 


1-D  FFT  method  with  n=200  @  t=1 


Figure  8.  1-D  FFT  method 
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c. 


COMSOL  MULTIPHYSICS  4.0A 


COMSOL  is  one  of  the  popular  computer  simulation  software  programs  used  to 
model  and  simply  translate  real-world  physical  laws  into  the  real  world  in  the  virtual 
form  using  the  finite  element  method.  COMSOL  is  a  commercial  problem-solving  tool 
that  produces  results  quickly.  However,  it  is  more  important  to  investigate  the  accuracy 
of  the  numerical  results.  We  compare  the  solution  differences  between  MATALB  and 
COMSOL  in  section  D. 

Figure  9  is  the  result  from  COMSOL  and  there  is  no  surprise  that  both  results 
from  MATLAB  and  COMSOL  are  very  close.  The  step-by-step  process  of  creating  this 
1-D  dithering  laser  problem  using  (Eq.  3.4)  in  COMSOL  is  attached  in  Appendix  C. 


COMSOL  1-D 


Figure  9.  1-D  COMSOL 
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D. 


COMPARISON  ON  A  MODEL  PROBLEM 


We  compare  the  results  from  the  1-D  Crank-Nicolson  method  and  COMSOL  and 
plot  the  relative  difference  in  Figure  10,  where  the  number  of  points  is  increased  from 
200  to  2001.  The  relative  error  between  MATLAB  and  COMSOL  results  is  about  10-4, 
which  is  really  small  and  tolerable.  Here,  we  only  show  the  comparison  between  Crank- 
Nicolson  method  in  MATLAB  and  COMSOL,  and  the  comparison  in  FFT  method  and 
COMSOL  are  also  pretty  similar.  After  several  tryouts,  it  is  intuitive  to  see  that  the 
relative  error  goes  down  once  we  increase  the  number  of  points  N,  but  it  takes  longer  to 
do  so. 


Differnce  between  C-N  Matlab  &  COMSOL  N=2001 


Figure  10.  Relative  error  plot  in  1-D  Crank-Nicolson  and  COMSOL 


E.  REAL  PROBLEM  SIMULATION  (STEEL  AISI  4340) 

After  we  have  demonstrated  that  both  MATLAB  and  COMSOL  are  giving  us 
acceptable  results,  we  can  choose  a  specific  material  to  simulate  a  real  problem.  The 
material  we  use  in  our  1-D  plot  comparison  is  Steel  AISI  4340,  a  built-in  material  in 
COMSOL  which  has  the  material  contents  listed  in  Table  2: 
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Property 

Name 

Value 

Unit 

Heat  Capacity 

475 

J/(kg*K) 

Density 

P 

7850 

Kg/(mA3) 

Thermal  Conductivity 

Kt 

44.5 

W/(m*K) 

Thermal  Diffusivity 

OCj 

Kt 

P'Cp 

mA2/s 

Melting  Point 

1783 

K 

Table  2.  Thermal  property  of  Steel  AISI 4340 


The  dithering  laser  we  deploy  has  the  inputs  in  Table  3: 


Laser  Input 

Name 

Value 

Unit 

Magnitude  of  Gaussian  source 

h 

1.0e9 

W/(mA3) 

Effective  radius  of  Gaussian  heat  source 

d 

0.02 

m 

Rotating  frequency 

Freq 

1 

Hz 

Laser  stop  time 

Time 

1 

S 

Center  of  rotation 

x0 

0.5 

m 

Rotating  radius 

a 

0.25 

m 

Table  3.  1-D  dithering  laser  input  on  Steel  AISI  4340 

Figure  11  from  MATLAB  and  Figure  12  from  COMSOL  depict  the  results  of 
temperature  rise  using  the  contents  from  Table  2  and  Table  3.  Both  results  are  very  close. 
Figure  12  is  from  the  Crank-Nicolson  method  and  the  FFT  method  yields  a  very  similar 
result.  It  is  clear  that  the  maximum  temperature  rise  decreases  as  the  rotating  period 
decreases;  in  other  words,  the  higher  the  frequency  is,  the  less  the  temperature  rise 
increases.  At  t  =  \s ,  Figure  13  shows  the  maximum  temperature  rise  of  the  steel  AISI 
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4340  versus  the  frequency  (reciprocal  of  the  period)  of  the  rotating  Gaussian  beam  from 
1  Hz  to  100  Hz  .  It  is  observed  that  a  higher  frequency  leads  to  lower  the  maximum 
temperature  rise.  It  should  be  pointed  out  that  our  results  are  consistent  with  earlier 
analytical  studies  for  the  semi-infinite  domain  [5]. 

Figure  14  shows  the  temperature  rise  as  a  function  of  time  at  a  fixed  point 
x  =  0.75  of  the  material  with  dithering  laser  beam  period=0.1  sec.  The  temperature 
increases  as  the  dithering  laser  beam  moves  closer  to  the  point  and  the  temperature 
almost  stays  steady  as  the  laser  beam  moves  away.  The  overall  line  shape  behaves  as  an 
increasing  function  of  time. 

Figure  15  shows  the  maximum  temperature  rise  versus  frequency  with  two 
different  heat  sources  70 .  The  higher  I0 ,  the  higher  the  temperature  rise. 


x 

Figure  1 1 .  Temperature  rise  on  the  material  of  steel  AISI  4340  using  the  1-D  Crank- 
Nicholson  method 
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I'D  CQMSOl 


Figure  12.  Temperature  rise  on  the  material  of  steel  AISI  4340  using  1-D  COMSOL 


Max  Temp  rise  VS  Frequency  at  time=1s 


Figure  13.  1-D  maximum  temperature  rise  of  steel  AISI  4340  versus  the  frequency  of 

the  dithering  laser  beam 
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Fixed  point  x=0.75  temperature  change 


Figure  14.  1-D  temperature  change  with  time  at  a  fixed  point  x=0.75  of  steel  AISI 

4340  with  dithering  laser  beam  period=0.  Is 


Max  Temp  rise  VS  Frequency  at  time=1s 


Figure  15.  The  maximum  temperature  rise  of  steel  AISI  4340  versus  frequency  with 
two  different  heat  source  I0 
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IV.  NUMERICAL  SOLUTION  FOR  A  TRANSIENT,  TWO- 
DIMENSIONAL  TEMPERATURE  DISTRIBUTION  IN  A  FINITE 
FILM  DUE  TO  A  ROTATING  OR  DITHERING  LASER  BEAM 


Now  we  extend  our  work  in  1-D  to  2-D. 


A.  THE  CRANK-NICOLSON  METHOD 


We  start  with  the  2-D  nonhomogeneous  heat  equation 


du 


=  aT 


82u  d 


dt 

f(x,y,t) 


-+- 


dxT  dy  ) 


a 


+-^f(x,y,t)  ,  0<x<Lx,  0 <y<Ly 
Kr 


Ind2 


xc(t)  =  x0  +a-cos 


yc(t)=y0+b-sin 


-yx-xyty+iy-yyty) 

,  2d 2 


»^=^T 


(Eq.  4.1) 


L 


(Eq.  4.1)  can  be  illustrated  as  in  Figure  16.  The  heat  source  f(x,y,t )  created  by  the  laser 
beam  is  illustrated  in  Figure  17. 


Oc(0,  y  ,(t)) 


(0,0)  (Lx, 0) 

Figure  16.  2-D  schematic  of  the  laser  beam  and  the  finite  work  piece 


35 


Figure  17.  Laser  beam  creates  heat  source  as  a  Gaussian  distribution  [After  2], 


We  discretize  this  2-D  rectangular  domain  as  shown  in  Figure  18.  In  order  to 
facilitate  our  presentation,  we  setZx  =Ly=  1  and  the  thermal  properties  aT  and  KT  equal 

du 

1,  thus  the  heat  equation  is  simplified  to  the  form  —  =  A u  +  f  ( x,y,t )  . 

dt 


Figure  18.  Discretization  of  a  2-D  rectangular  region 
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The  discretized  point  is  denoted  as  (x. ,  yj )  ,  where  i  goes  from  1  to  n  and 
j  varies  from  1  to  m  .  Here,  n  and  m  are  the  total  number  of  grid  points  in 
x  and  y  respectively.  More  explicitly, 


x.  = 


f.  o 

l  — 

V  2  j 


Ax. 


,  Ax  =  -^ - -,  i  =  l,2,...n 


y,= 


J- T 

V  2  j 


Af 


(Ly-  0) 

,  Ay  =  - - j  =  l,2,...m 

m 


The  boundary  conditions  assume  that  the  solid  is  insulated  at  the  edges 
(boundaries)  and  the  initial  condition  is  that  there  is  no  temperature  rise  with  respect  to 
the  ambient  temperature  initially.  Thus,  we  have 


BCs: 


du 

dx 

du 

¥ 


x=0 


y= o 


du 

dx 

du 

dy 


X=1 


y= i 


=  0 

=  0 


ICs:  u(x,y,  0)  =  0 


We  introduce  the  shorthand  notation  uki  j  =  u(xi,yj,tk) .  Then  the  Crank-Nicolson 
2-D  scheme  of  (Eq.  4.1)  becomes 


k+ 1  .  k 


U..  -«  1 


At 


+ 


+ 


Ui+\J  ~^UiJ  +Ui-l,j  ,  C  -24+1+<t 


Ax2 


+- 


Ax2 


,/  -Pm*  +M*  >/+1  -2w*+1+z/*+1 

”  j+i  ZMi j  +  Mi j-i  ,  ”  j+i  j  +  j-i 


-  + 


A/  Ay2 

^  [/(*,  ,yj,tk)+ f(xt ,  y} ,  tk+l )] 


(Eq.  4.2) 


1  <  i  <  n,  1  <  j  <  m 
k  =  0,1,2....  tk=kAt 


The  BC  can  be  satisfied  by  setting 
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k 

k 

UKi 

=  u0  J 

k 

k 

u  . 

n,j 

=  Un+lJ 

k 

k 

uli 

=  Ui,  0 

k 

_  k 

u. 

i,m 

~  Ui,m+\ 

We  can  rearrange  the  temperature  at  the  interior  grid  points  to  form  a  vector  u  as 
depicted  in  Figure  19: 


Figure  19.  Putting  the  temperature  of  each  point  in  a  vector  form  u  after  discretization 


We  will  convert  the  linear  system  (Eq.  4.2)  into  the  following  matrix  equation: 

-*+ 1  —  - 

Au  =  Bu  +  f=  b 
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where  A  and  B  are  the  matrices  to  be  found,  [  u  I  are  unknown  and  j  w  1  are  known. 


We  rewrite  (Eq.  4.2)  as 
k  At 

Ui,  =u\ ,  + - - 

,J  hJ  2  Ax2 


[umj  ~  2uu  +  Ui-lJ ]  +WVj[uiJ+ 1  _  2ul  +  UiJ~l ] 


2  Ay 


.  ^  r  k+i  o  it+I  .  fc+i  "I  .  At  r  £+1  0  £+1  .  £+1  1 

+  “  “  [UMJ  -  2ui,j  +  ui-u  J  +  L“«v+1  _  2“i.y  +  Mw-1  J 


+ 


2Ax2 
At 


—  [/(A- >  A,'  v 4  )  +  ./'(A ,  JO ,  4+1 )] 


Putting  m  on  the  left-hand  side  of  the  equation,  a  and  /  on  the  right-hand  side, 
we  obtain: 


At 


u^'  + 


2 Ax2  MJ 

At  /t 

"W,. . ,  + 


2Ax 

At 


2  i+l, y 


,  At  At 

i+ — r+ — 7 

Ax2  Ay2 

At  At 


k+l_  &  kl 
l’J  t  a  „2  i-l,y 


1- 


Ax2  Ay2 


2Ax 
At 


2 Ay2  ,J+1  2A y2 


At 


At 


W,.  ■  H - —  W .  i  j  H - —  W/  H - -  Uj 


2 Ax2  ’~u  2 Ay2  /J+1  2Ay2 


(Eq.  4.3) 


+ — [/(a  > jo  > 4 )  +  /(a  > jo >  4+i )] 


We  start  a  simple  case  with  i  =  1,  j  =  1  first, 


At 

-w, ,  + 


2Ax 


2  2,1 


,  At  At 

1  +  — r  + 


At 


2Ax 


At, 


Wt,  + 


2  2,1 


1- 


Ax2  Ay2 

At  At 


n,i 


At 


2Ax 


,k+l 

2  U0,l 


At 


£+1 


At 


,-W,  2 

2Ay2  1,2  2A_y 


2  Wl,0 

,fc+l 


Ax2  Ay2 


At 


At 


At 


m(i  H - Twni  "I - irwf,  "I - Tu\(\ 

1,1  2Ax2^  2 Ay2  1,2  2Ay2  ^ 

=M1,1 


=«1,1 


+  y[/(A>Ai,4)  +  /(A > A, ,  4+1 )] 


Using  the  boundary  conditions  (  =  m**1  ,  and  so  on),  we  can  simplify 

this  equation 
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,  At  At 

1  + - r  +  - 


2  Ax  2Ay2 

At  At 


k+l 


At 


k+\ 


At 


u  2Ax2  24  2Ay 


ik+l 

2  Ul,2 


2Ax2  2A y2 


At  ^  At  fc  At  v  r /  \  r/  \~\ 

M,l+2VMl’2 


U\,\  ^  A  2  w2,l 


2Ax 


— t+i  ->*  — 

Recall  the  matrix  form  riu  =Bu  +  f  .  We  can  obtain  the  corresponding 
components  of  ri  and  5  from  this  equation: 


_  IX  1  At  At 

^(M)  -  1+  _  7  2  + 


-4(1,2)  =  - 


2 Ax2  2Ay2 
At 


2Ax2 


^(1,  »  +  !)  =  - 


At 


5(1,1)  =  1- 


2Ay2 

At  At 


5(1,2)  = 


2  Ax  2Ay2 

At 


2Ax2 


B{\,n  +  \)  = 


A  t 

2V 


and  we  will  put  y  [/(•*1,.y1,tit)-l-  /(x^y,,^,)]  into  a  component  of  a  vector  function 


called  /. 


When  i  =  2,  j  =  \,  (Eq.  4.3)  becomes 


At 


2Ax 


u*!1  + 


2  3,1 


,  At  At 

1  + - 7  + - r- 

Ax  Ay 


A t  ,,/t+i  At  ,,/t+i  At  ,,/t+i 

i  9  ^9  2  9  ^2  0 

2A_y2  2,2  2Ay2  ^ 

=U7  1 


A:+l 

24  2 Ax2  “u 


At 


2Ax 


At, 


wf,  + 


2  3,1 


At  At 


Ax2  Ay2 


^  At  ^  At  ^  At  ^ 

2,1  2 Ax2  1,1  2 Ay2  2’2  2Ay2 

=W2,1 


+  y  [/(*2 > 7i  v 4 )  +  /(*2 > Ti > 4+! )] 


To  simplify,  this  is  equivalent  to 
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At  fc+i 
-u\f  + 


2  "3,1 


2Ax 

~M2, 3) 


At  k 
u a  ,  + 


1  H - y  H - y 

Ax2  2Ay2 

^02) 
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2AT 
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2V 
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2Ax*  3,1 

5(2,3) 


A^  A^ 


Ax2  2Ay2 

5(2,2) 
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At 

2A?  "u 

5(2,1) 


A^ 


2Ay 

5(2,n+2) 


2  W2,2 


+ y  [/( Wi ,  4 )  + /(*2 > J4 > 4+i )] 


When  i  =  3,  j  =  1 ,  (Eq.  4.3)  becomes 


A/ 


2Ax 

A? 


w!!1  + 


2  4,1 


,  A/  A/ 
1  +  — ^  + 


2Ax 


A?, 
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2  4,1 
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Ax2  Ay2 
At  At 


,.k+\ 

*3,1 


At 
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At 


- u 0 


X+l 


A^ 


Ax2  Ay2 


Wg  1  + 


2 Ax2  2,1 

k 


2  Ay2” 1/3,2  2Ay 


2  ^3,0 

— m3,1 


A^  ^  At  ^ 

y  1  H  y  j  “I  y  a 

2 Ax2  2,1  2 Ay2  3’2  2Ay2  X2 

=w3,l 


+  y[/(x3  >yi>  4 )  +  /(x3 > v, ,  4+1 )] 


To  simplify,  this  is  equivalent  to 


At 


2Ax 

yy 


hM1  + 


2  4,1 


,  A^  A^ 

1  H - y  H - y 

Ax2  2Ay2 
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2Ax^ 
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2  4,1 


1- 


-4(3,3) 

At  At 


2Ax 

yyr 


2Ay 

-4(3,  «+3) 


/+1 
2  "3,2 


Ax2  2  Ay2 


£  A^  £  A^  At  r  y>/  \  rs  \~\ 

1/3,1  +  2A? 1/2,1  +  2^^ 1/3,2  +^^X35^15^+^X35’>;i’^+1^ 


5(3,3) 


5(3,2) 


5(3, n+3) 


When  i  =  n,  j  =  1 ,  (Eq.  4.3)  becomes 
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2Ax 
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Ax2  Ay2 

At  At 
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At 
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k+ 1 


2 Ax2  n‘u  2Ay2  n’2  2Ay2  ^ 


Ax2  Ay2 


£  A^  u  At  k  At 
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+  y  [/(*„ > Tp  4 )  +  f(xn > Ti ,  4+1 )] 
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To  simplify,  this  is  equivalent  to 


,  At  At 

1  + - r  +  - 


2Ax 2  2A y2 
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At 
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At 
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Now  i  =  1,  j  =  2 ,  (Eq.  4.3)  becomes 
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At 
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+  y  [/(*i ,  T2 >  ■ 4  )  +  /(*i > T2  >  4+i )] 


To  simplify,  this  is  equivalent  to 
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The  same  arguments  hold  for  i  =  1  and  j  =  3 . 
For  i  =  1,  j  =  m ,  (Eq.  4.3)  gives 
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To  simplify,  this  becomes 
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For  i  =  2,  j  =  m , ,  (Eq.  4.3)  takes  the  form 
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To  simplify,  this  becomes 
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Finally,  when  i  =  n,  j  =  m  ,  (Eq.  4.3)  becomes 


43 


At 

-ukl + 


2  Ax 


2  n+\,m 


,  At  At 

1  +  — 7  +  — 7 

Ax2  Ay2 


ut  t!  - 


At 


k+ 1 


At 


k+ 1 


At 


k+ 1 


/7,m  ^  *  2  n-\,m  ^  A  2  «,m+l  ^  *  2  «,m-l 

2 Ax  2 Ay  ■ — * — ■  2Ay> 


A? 


2Ax 


„  + 


2  «+l,m 


1- 


At  A? 


Ax2  A/ 


^  At  ^  At  ^  At  ^ 
w„,m  +T7~rl4-i,m  +rrr"»,»ti  +^rrrw«,m-i 


2Ax 


2A/ 


2A/ 


+ y  [/(*„ ’ ,  4 ) + /  (a,  ,  Tm ,  4+i )] 


To  simplify,  this  becomes 
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We  can  put  all  the  components  of  4  and  5  together  to  obtain  the  two 
matrices  4(„.m>n.m)  and  5(„.m,n.m) . 

-k+\  -k  — 

After  obtaining iwi  =  Bnmxnmunm,i  +  f  nmxl  =bnm,i,  we  solve  it  in  MATLAB 

-~k+ 1 

by  w  =  A\b  . 

We  use  (Eq.  4.1)  and  the  parameters  in  Table  4  to  show  the  temperature  rise  of  a 
model  problem  with  the  numerical  points  n  =  m  =  128  in  Figure  20.  The  MATLAB  code 
is  attached  in  Appendix  D. 


’4+i)] 
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Figure  20.  Snapshots  of  the  temperature  rise  on  a  film  of  a  model  problem  induced  by 
a  rotating  Gaussian  beam  using  the  2-D  Crank-Nicolson  method 


B.  THE  FAST  FOURIER  TRANSFORM  (FFT)  METHOD 

The  structures  of  Fast  Fourier  Transform  in  1-D  and  2-D  are  similar;  instead  of 
using  fft  and  inverse  Fast  Fourier  Transform  ifft  commands  in  1-D  MATALB  code,  the 
2-D  code  uses  fft2  and  ifft2  to  carry  out  the  computation.  Briefly  speaking,  one  needs  to 
even  extend  the  problem  from  domain  [0,l]x[0,l]  to  [0,2]x[0,2],  make  the  problem 
periodic  in  both  x  and  y  direction  with  period  2,  and  then  apply  the  Fourier  transform  and 
its  inverse  to  obtain  the  numerical  solution.  Our  MATLAB  code  is  attached  in  Appendix 
E. 
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Figure  21  uses  the  input  parameters  from  Table  4  and  returns  a  very  similar  result 
to  Figure  20.  The  2D  FFT  has  a  better  performance  than  the  2-D  Crank-Nicolson  method. 
In  other  words,  FFT  can  produce  the  result  faster  than  the  Crank-Nicolson  method  with 
the  same  number  of  numerical  points.  We  will  discuss  the  details  in  Section  D. 


(Figure  continued  on  next  page.) 
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Figure  2 1 .  Snapshots  of  the  temperature  rise  on  a  film  of  a  model  problem  induced  by 
a  rotating  Gaussian  beam  using  2-D  FFT  method  in  3-D  view 


C.  COMSOL 

Figure  22  (a)  is  the  result  from  COMSOL  and  Figure  22  (b)  is  the  result  using  the 
Crank-Nicolson  method  in  MATALB  based  on  (Eq.  4.1)  and  Table  4  with  the  numerical 
points  n  =  m  =  128  in  x  and  y  direction,  respectively.  The  results  are  very  close  and  we 
will  do  a  point-by-point  error  analysis  in  Section  D.  The  detailed  process  of  creating  this 
2-D  COMOSL  is  attached  in  Appendix  F. 
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Figure  22.  Model  problem  comparison:  (a)  COMSOL  (b)  MATLAB  FFT  2  method 


D.  COMPARISON  ON  A  MODEL  PROBLEM 

As  we  have  mentioned  in  Section  B,  FFT  returns  the  solution  faster  than  Crank- 
Nicolson;  the  efficiency  comparison  is  shown  in  Figure  23.  When  N=29  points,  we  can 
see  that  FFT  takes  about  10  seconds  but  Crank-Nicolson  takes  more  than  10  seconds  to 
finish  the  computation.  So  FFT  is  about  ten  times  faster  than  the  Crank-Nicolson  method 
for  this  test  problem. 
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Figure  23.  Efficiency  plot:  Crank-Nicolson  2-D  method  versus  FFT  2-D  method 


In  Figure  24,  by  using  the  result  from  Figure  22,  we  compare  the  relative 
difference  in  temperature  rise  with  respect  to  time  at  a  fixed  point  (x=0.5,  y=0.5)  in 
MATALB  FFT  2  method  and  COMSOL  using  number  of  point  N=256.  The  relative 

'J 

error  is  about  10'  and  this  result  is  tolerable.  The  relative  error  goes  down  as  we  increase 
the  number  of  numerical  grid  points  N. 


Figure  24.  Relative  error  plot  in  2-D  FFT  method  and  COMSOL 
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E.  REAL  PROBLEM  SIMULATION  (STEEL  AISI  4340) 


We  use  Steel  AISI  4340  as  our  test  material  with  material  properties  from  Table  2 
and  the  rotating  laser  we  deploy  has  the  following  input  in  Table  5: 


Input 

Value 

Unit 

h 

1.0e8 

W/mA3 

d 

0.02 

m 

4 

1 

m 

4 

1 

m 

*0 

0.5 

m 

To 

0.5 

m 

a 

0.25 

m 

b 

0.25 

m 

period 

1 

s 

Table  5.  2-D  rotating  laser  input  on  Steel  AISI  4340. 

Figure  25  depicts  the  temperature  rise  at  different  times  within  one  period.  It  is 
observed  that  heat  will  not  spread  out  quickly  enough  due  to  the  properties  of  the  Steel 
AISI  4340  material  so  the  temperature  at  those  points  directly  shined  by  the  laser  rise 
quickly.  However,  those  points  far  away  from  the  points  hit  by  laser,  for  instance,  the 
center  point  (x=0.5,  y=0.5),  has  almost  no  temperature  rise. 
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Figure  25.  Snapshots  of  the  temperature  rise  on  a  Steel  AISI 4340  film  induced  by  a 
rotating  Gaussian  beam  using  COMOSL  with  period=ls 


Figure  26  shows  the  temperature  rise  as  a  function  of  time  at  a  fixed  point 
(x=0.75,  y=0.5)  of  Steel  AISI  4340  with  dithering  laser  beam  period=0.1  sec.  The 
temperature  increases  as  the  laser  beam  rotates  close  to  the  point  and  the  temperature 
almost  stays  steady  as  the  laser  beam  moves  away.  The  overall  line  shape  behaves  as  an 
increasing  function  of  time,  as  expected. 
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Fixed  point  (x=0 .75  y=0.5)  temperature  change 


Figure  26.  2-D  temperature  change  with  time  at  a  fixed  point  of  Steel  AISI 4340  with 

rotating  beam  period=0.1s 


The  quantitative  relationship  between  the  maximum  temperature  rise  and  the 
rotating  frequency  at  time=ls  is  depicted  in  Figure  27;  the  maximum  temperature  rise  is  a 
decreasing  function  of  the  frequency  (reciprocal  of  the  rotating  period). 


Max  Temp  rise  VS  Frequency  at  time=1s 


Figure  27.  2-D  maximum  temperature  rise  of  Steel  AISI  4340  versus  the  frequency  of 

the  rotating  laser  beam 


53 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


54 


V.  NUMERICAL  SOLUTION  FOR  A  TRANSIENT,  THREE- 
DIMENSIONAL  TEMPERATURE  DISTRIBUTION  IN  A  FINITE 
SOLID  DUE  TO  A  ROTATING  OR  DITHERING  LASER  BEAM 


We  have  already  verified  that  COMSOL  has  returned  a  very  accurate  solution 
compared  with  other  numerical  methods  in  both  1-D  and  2-D  codes.  Therefore,  instead 
of  writing  a  complicated  MATLAB  3-D  code,  we  use  COMSOL  to  obtain  the  3-D 
answer.  Recall  (Eq.  4.1)  and  impose  insulated  boundary  conditions.  We  use  Steel  AISI 
4340  as  our  test  material  with  material  properties  from  Table  2  and  the  rotating  laser  we 
deploy  has  the  following  input  from  Table  6: 


Input 

Value 

Unit 

h 

5.0e5 

W/mA2 

d 

0.02 

m 

4 

1 

m 

4 

1 

m 

4 

1 

m 

*0 

0.5 

m 

To 

0.5 

m 

a 

0.25 

m 

b 

0.25 

m 

period 

1 

s 

Table  6.  3-D  rotating  laser  input  on  Steel  AISI  4340 

Figure  28  depicts  the  temperature  rise  at  different  times  within  one  period.  The 
hottest  spot  is  where  the  point  hit  instantly  by  the  laser  as  in  the  1-D  and  2-D  cases.  The 
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maximum  temperature  rise  is  1468K.  Those  points  far  away  from  the  area  hit  by  the 
laser,  for  instance,  the  center  point  (x=0.5,  y=0.5,  z=l),  have  little  temperature  rise. 


Figure  28.  Snapshots  of  the  temperature  rise  on  a  Steel  AISI  4340  solid  induced  by  a 
rotating  Gaussian  beam  using  COMSOL  with  period=ls 


Figure  29  shows  the  temperature  rise  at  a  fixed  point  (x=0.75,  y=0.5,  z=l)  as  a 
function  of  time.  Figure  30  shows  the  maximum  temperature  rise  of  the  whole  domain  as 
a  function  of  time;  the  overall  hottest  spot  is  at  (x=0.341,  y=0.307,  z=l)  when  time=0.64s 
with  the  rotating  laser  beam  period=ls.  Figure  31  and  Figure  32  depict  the  temperature 
rise  at  fixed  time=ls  in  different  layers.  It  is  observed  that  the  heat  does  not  spread 
downward  quickly  and  there  is  almost  no  temperature  rise  0.1m  below  the  top  surface 
shined  by  the  laser.  After  several  tryouts,  materials  have  larger  values  in  diffusivity  than 
Steel  AISI  4340  can  make  heat  spread  out  faster  and  so  make  the  temperature  rise  higher 
than  Steel  AISI  4340. 


'-,U' 


--,U 


56 


□  D.2  o.4  o.e  o.e  i 

fimsfe) 


Figure  29.  3-D  temperature  change  as  a  function  of  time  at  a  fixed  point  of  steel  AISI 

4340  with  rotating  laser  beam  period=ls 


Max  temp  with  rotating  laser  beam  period- Is 


Figure  30.  Maximum  temperature  as  a  function  of  time  with  rotating  laser  beam  period 
=1  s 
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Figure  31.  3-D  temperature  rise  at  different  layers  from  z=0  to  z=l 
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Figure  32.  Maximum  temperature  at  different  depth  (a)z=l  top  surface  (b)z=0.99 
(c)z=0.95  (d)z=0.90 
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Figure  33  shows  the  maximum  temperature  rise  at  time=ls  when  the  period  is 
reduced  to  0.1s;  the  laser  beam  rotates  10  cycles.  The  maximum  temperature  rise 
decreases  from  1754K  of  period=ls  to  670K  of  period=0.1s.  Therefore,  the  maximum 
temperature  rise  can  be  reduced  by  increasing  the  frequency  of  the  rotating  laser  beam. 


Figure  34  shows  the  temperature  rise  at  a  fixed  point  (x=0.75,  y=0.5,  z=l)  of 
period=0.1s  as  a  function  of  time.  The  overall  temperature  rise  oscillates,  and  its 
envelope  behaves  as  an  increasing  function  of  time,  but  its  maximum  temperature  rise  is 
smaller  compared  to  the  maximum  temperature  rise  with  period=ls,  as  illustrated  in 
Figure  29.  Figure  35  shows  the  maximum  temperature  rise  of  the  whole  domain  as  a 
function  of  time,  and  its  envelope  behaves  as  an  increasing  function  as  well.  Figure  36 
depicts  the  quantitative  relationship  between  the  maximum  temperature  rise  and  the 
rotating  frequency  at  time=ls.  The  results  agree  with  earlier  analytical  studies  for  the 
semi-infinite  domain  [5]. 


time -0.1s.,  perlotl-'O.ls 
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Figure  33.  Snapshots  of  the  temperature  rise  on  a  Steel  AISI  4340  solid  induced  by  a 
rotating  Gaussian  beam  using  COMSOL  with  period=0.1s 
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Fi*ed  point  <k=Q.75.  y-o  5.  z=-i>  change 


Figure  34.  3-D  temperature  change  as  a  function  of  time  at  a  fixed  point  of  Steel  AISI 

4340  with  rotating  laser  beam  period=0. 1 


Max  tamp  with  rotating  bean  Deriod=0.  Is 


Figure  35.  Maximum  temperature  as  a  function  of  time  with  rotating  laser  beam  period 
=0.1s 
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Figure  36.  3-D  maximum  temperature  rise  of  Steel  AISI 4340  versus  the  frequency  of 

the  rotating  laser  beam 
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VI.  CONCLUSIONS  AND  FUTURE  WORK 


In  this  thesis,  both  analytical  and  numerical  solutions  for  describing  the  transient 
temperature  rise  induced  by  a  moving  laser  in  a  finite  domain  have  been  developed.  We 
have  exploited  several  methods  including  the  eigenfunction  expansion,  the  Crank- 
Nicolson  scheme,  FFT  and  COMSOL. 

We  have  confirmed  that  the  faster  the  laser  rotates  (i.e.,  the  higher  the  frequency) 
the  lower  the  temperature  rise  induced.  In  other  words,  to  reduce  the  military’s 
vulnerability  to  high-energy  laser  weapons  it  is  possible  to  let  the  object  rotate  or  rock  to 
minimize  the  temperature  rise.  The  quantitative  relationship  between  maximum 
temperature  rise  and  laser  beam  rotation  frequency  can  be  used  as  a  general  guide  for 
adjusting  the  speed  of  rotation  of  the  object  in  order  to  prevent  temperature  rise  from 
reaching  the  melting  point. 

This  thesis  can  be  explored  deeper  in  the  future.  Some  future  potentials 
endeavors  include  but  not  limited  to: 

A.  Increase  or  decrease  the  effective  radius  of  the  laser  beam  d  in  Figure  3 
and  Figure  17  to  analyze  how  temperature  rise  is  affected. 

B.  Increase  or  decrease  the  radii  a  and  b  of  the  rotating  trajectory  of  the 
Gaussian  beam  in  Tables  3,  5  and  6  to  analyze  how  maximum  temperature 
rise  is  affected. 

C.  Create  3-D  MATALB  codes  and  compare  the  results  with  COMSOL. 

D.  Instead  of  using  TEM^  mode  Gaussian  distribution  as  the  heat  source 
illustrated  in  Figure  17,  different  transversal  modes  in  a  laser  spot  such  as 
TEM0l  or  TEMn  ,  can  be  used  to  further  seek  the  analytical  and  numerical 
solutions  [8], 

E.  Thermal  properties  are  assumed  to  be  temperature  dependent:  this  makes 
the  nonhomogeneous  heat  equation  nonlinear,  but  more  realistic  [3].  This 


63 


nonlinear  model  can  be  solved  without  mathematical  background  in 
nonlinear  programming  using  appendix  G,  COMSOL  code  for  3-D 
simulation  by  imposing  certain  materials  whose  thermal  properties  are 
temperature-dependent,  such  as  silver  [9], 

F.  In  COMSOL  3-D  geometry,  try  cylindrical  coordinate  and  spherical 
coordinate  rather  than  Cartesian  coordinate. 

G.  Conduct  an  experiment  and  see  if  the  theoretical  modeling  is  accurate. 
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APPENDIX  A.  CRANK-NICOLSON  CODE  FOR  1-D  SIMULATION 


N=256; 

%  number  of  numerical  points  in  [ 

dx=l/N; 

%  spatial  step 

x— ( [ 1 : N] ’-0.5) *dx; 

%  numerical  grid 

u=zeros (N, 1 ) ; 

%  solution  at  current  time,  u(j): 

d=  0 . 0  2  ; 

%  radius  of  Gaussian  source 

a=l . 0 ; 

%  magnitude  of  Gaussian  source 

u0=zeros (N, 1) ; 

%  initial  value 

u=u0 ; 

o 

o 

dt=dx/ 8 ; 

%  time  step 

T=2  5/ 256; 

m=T/dt ; 

%  number  of  time  steps 

4-> 

rs 

-x 

e1 

o 

II 

^3 

o 

o 

r=dt/ dxA2  ; 

%  forming  the  matrices 

d0= [ 0 . 5 ;  ones (N-2 ,  1 )  ;  0.5]*r; 

dl=0 . 5* ones (N, 1 ) *r  ; 

d_l  =  0 . 5* ones (N, 1 ) *r; 

A=spdiags ( [ -d_l ,  1+dO,  -dl ]  ,  [-1,0,1] ,N,N) ; 

B=spdiags ( [d_l ,  1-dO,  dl ] ,  [ -1 , 0 , 1 ] ,  N,  N)  ; 

o 

o 

y0  =  0 . 5+0 . 25*sin ( 10*2*pi*0 ) ;  %  location  of  Gaussian  source  at  tA{k-l} 

f 0=a*exp ( - (x-yO ) . A2 / ( 2 *dA2 ) ) /sqrt ( 2 *pi*dA2 ) ;  %  Gaussian  source  at  tA{k- 

1 } 

o 

o 

plot (x, u,  1 b- ' ,  1 linewidth  *,2.0) 
axis ( [0,1,-0.04,0.16]  ) 
drawnow 
for  k=l : m, 

yl  =  0 . 5  +  0 . 25*sin  (10*2*pi*t  (k+1)  )  ;  %  location  of  Gaussian  source  at  tAk 

f l=a*exp (- (x-yl) . A2/ (2*dA2) ) /sqrt (2*pi*dA2) ;  %  Gaussian  source  at  tAk 

b=B*u+dt* ( f 0+f 1 ) /2 ;  %  right-hand-side  of  the  linear  eq 

for  uA{k} 

u=A\b;  %  solving  the  linear  eq  for  uA{k} 

f  0 = f  1 ; 

%pause (0.1) 

plot(x,u, 1 b- ' , 1 linewidth ',2.0) 
axis ( [0,1,-0.04,0.16]  ) 
drawnow 
end 
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APPENDIX  B.  FFT  CODE  FOR  1-D  SIMULATION 


N=256; 

dx=l/N; 

x— [ 0 : N ]  ' *  dx  ; 

u=zeros (N+l, 1) ; 

d=0 . 02; 

a=l; 

u0=zeros (N+l, 1) ; 
u=u0 ; 

o 

o 


%  number  of  numerical  points  in  [0,  1] 

%  spatial  step 
%  numerical  grid 

%  solution  at  current  time,  u(j)=u(x(j)) 
%  radius  of  Gaussian  source 
%  magnitude  of  Gaussian  source 
%  initial  value 


dt=dx/8;  %  time  step 

T=25/ 256; 

m=T/dt;  %  number  of  time  steps 

t= [ 0 : m] *dt ; 

o 

o 

%  going  to  the  coefficients  of  cosin  expansion 
w= [ u ;  u ( N : - 1 : 2 ) ] ; 
z=fft (w) ; 

cu=real (z (1 :N+1) ) /N; 
r= ( [0 :N] 1 *pi) . A2 ; 

hl= [dt ;  (1-exp (-r (2 :N+1) *dt) ) . /r (2: N+l) ] ; 

h2= [ 0 . 5*dt A2 ;  ( 1- (r (2 : N+l ) *dt+l ) . *exp ( -r (2 : N+l ) *dt ) ) ./r(2:N+l) . A2] ; 

o 

o 


y0  =  0 . 5+0 . 25*sin ( 10*2*pi*0 ) ;  %  location  of  Gaussian  source  at  tA{k-l} 

f 0=a*exp ( - (x-yO ) . A2 / ( 2 *dA2 ) ) /sqrt ( 2 *pi*dA2 ) ;  %  Gaussian  source  at  tA{k- 

1 } 

w= [ f  0 ;  fO (N:-l:2)  ] ; 
z=fft (w) ; 

cf0=real (z (1 :N+1) ) /N; 

o 

o 


plot(x,u, ' r- 1 , ’ linewidth ',2.0) 
axis ( [0,1,-0.04,0.16] ) 
drawnow 
for  k=l : m, 

yl  =  0 . 5  +  0 . 25^sin  (10^2,lrpi^t  (k+1)  )  ;  %  location  of  Gaussian  source  at  tAk 

f l=a*exp ( - (x-yl ) . A2 / ( 2 *dA2 ) ) /sqrt ( 2 *pi*dA2 ) ;  %  Gaussian  source  at  tAk 

w= [ f 1 ;  fl (N : -1 : 2 ) ] ; 
z=fft (w) ; 

cf l=real ( z ( 1 : N+l ) ) /N; 

%  update  cu 

cu=cu .  ^exp  (-r*dt)  +cfl.*hl+  (cfO-cfl)  /dt.,lrh2; 
cf 0=cf 1 ; 

%  going  back  to  the  function 
z=N*[cu;  cu(N:-l:2)]; 
w=if f t ( z ) ; 
unreal (w ( 1 : N+l ) ) ; 

%pause (0.1) 

plot(x,u, ’ r- ' , 1 linewidth ',2.0) 
axis ( [0,1,-0.04,0.16]  ) 
drawnow 
end 
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APPENDIX  C.  COMSOL  CODE  FOR  1-D  SIMULATION 


1 .  Open  COMSOL  4.0a  with  ID  and  hit 


i  Model  Wizard  H  Model  Library 

°  □ 

Select  Space  Dimension 

0  ^ 

0  3D 

2D  axi  symmetric 

02D 

ID  axi  symmetric 
•  ID 
OD 

2.  In  Heat  Transfer,  select  Heat  Transfer  in  Solids  (ht)  and  right-click  Add  selected  then 
rit  . 

Add  Physics 

O  Recently  Used 

t>  AC/DC 
|>  i=3 '■  Acoustics 

!,  Chemical  Species  Transport 

O  ij _ ff  Electrochemistry 

l>  Fluid  Flow 
^  \S\  Heat  Transfer 

^  Heat  Transfer  in  Solids  (ht) 
ggji  Heat  Transfer  in  Fluids  (f  ^  Add  Selected 

c3jjl  Heat  Transfer  in  Porous  - 

s  Heat  Transfer  with  Surface-to-Surface  R.adiation  (ht] 

3  i  o  h  eat  Tra  n  sf  er  (ht] 

#  Surface-to-Surface  R.adiation  (rad] 

Joule  Heating  (jh] 
f:  Q  P I  a  sm  a 
:  Au  Mathematics 


^  X 

Selected  physics 

^  Heat  Transfer  (ht] 

3.  In  Preset  Studies,  select  Time  Dependent  and  hit  Finish  ^  . 

J  Model  Wizard  |®]  Model  Library:  ^  E3 

Select  Study  Type  <=>  i=>  |^~ | 


Studies 


Selected  physics 

±  H  eat  Tra  n  sf  er  (ht) 
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4.  Under  Model  Builder,  right-click  on  Global  Definitions  and  left-click  to  select 
Parameters,  input  Parameters  as  following: 


Mcde)  Buiict' 


Settings 


Pi  Parameters 

•  Parameter* 


*  ’  i  l-P.mphi'wfl 
t  =  Global  CViniticns 
Pi  Psremettrs 
4  ■  Model  lltoodJl 
p  =  Definit  jns 
!•  f\  G«m«ry  1 
£  Materials 
\  Hsat  T'ensterpiC 
P@M«hl 

*  r-  Study  1 

>Stept:frorr(Jtoperi5dit: 

;  |Sw  Sctotr  Configurations 

m  1  ‘i  RernNr 

5.  Right-click  on  Definitions  and  left  click  to  select  Variables,  input  Variables  as 
following: 


Name 

Expression 

Vjljt 

Description 

d 

0J2M 

OMm 

Gaussian  km  Radius 

aericc 

Its] 

Is 

i 

0.2S;m| 

0.23  m 

rangt  of  rotation 

>0 

D.S(ir| 

0.5 

center  of  rotation 

B 

lKW/m*3t 

lW/m 1 

magnitude  of  Gaussian  3eam 

del  Builder 


1-D.mph  vw£) 

S  Globe  I  Definitions 
Model  1  (modi! 

=  Definitions 
i "  Variables  1 
Q^Viwl 
\  Geoimetiyl 
$  Materials 
|  Hea:Transfer|'btj 
%  Mesh  I 

^  5h4'l 

6. 


Settings  .![  Model  Libraey 


Variables 

Geometric  Scope 
Geometric  entity  I  evt: 

^  Variables 


Entire  model 


Name 

Expression 

Unit 

2C 

Jd0+a*sin{2*pi  "t^pefiod] 

m 

f 

Wsqrt(2Tpi*d  ■r-2)[Lm]Katp(-(3i-H)  A2/g\P2)] 

W/m: 

Right-click  on  Geometry  and  left-click  to  add  Interval.  The  Left  endpoint  is  0  and 


the  Right  endpoint  is  1 .  Left-  click  the  build  all  button 
f  teote\  Bu  Idsr  v  °  □ 


4  *  1-Dmph  [jneofl 

B  Global  Dtfiiiittens 
4  Model  1  fmotfy 
B  Definitions 
a  A,  Geometry  I 
A  Interval  1  (il) 
Form  Union  [fir 
a  ©  Materials 


Settings  Model  Library 

Interval 

*  Interval 


to  create  geometry. 

iSI8Lttc 


Number  off  intervals;  One 
Left  endpoint: 

Right  endpoint: 


1 


7.  Right-click  on  Materials  and  left-click  on  add  Materials.  Make  all  the  values  in 
Material  Contents  equal  1.  Right-click  on  the  “line”  in  Graphics  and  left-click  on  it  to 
select.  Make  sure  that  “  1  ”  is  under  the  selection. 
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i  tic-*  *l  -;i  a  i..  tL  i  s  O*  ■ " 


■X4  Qd  -  ■  ■ 

El  **- 

a^.'«E 

Vfacri  I  'mil 
f-  £di  H  i>i 
}?■.  h;  nn  ■-.  I 

btniM  j'd, 

$  Fj  p«  LNw-fr  j 

*M|#Wllfc 

|  HWTyilk1  H1 

■  jig.1 .  I'wr  J  hjj 

V  l***!  TrrlijMV 
|  HlUll 
',  vdi’»b 
jr;~H.te  ttij« 

■J  Ti:  ■ 

U  ]DF1«i>llT  L 


t  .  rT  :  £r.::il.imv 


W  -  □  ct  -npp  n 


%  Material 

(m  iHvir  rnyf>k 


jccmlnki-H?  r  ri1  L'.'tbi 

*1 

"jIh-utv  K'v.  * 

__  '\ 

▼  '.’dlcnd  ^wrthi 


■  >it  rM-. 

C4  FW  I 


E-rr:  fizirma 
BtTKtwttjy 
[  hncTuq  vi  r  Pi' = J  * 

S-i :  MkaJhrii :: 
FyjMfcdSf  W?1W 
uiiHiiu 


■■  VdllTj  IiTfi  Mi 


Jr^cTi  ?ikm 

■Jr.  j  ■'Tjwt.  aa.n 

TtafrVltMtfcrfvlIV 

x  L 

'AGO  l»k 

Iteisty 

tx  1 

tj’V|  |W 

Par  l  riirl  :iiiji 

<4  i 

1  fc=“F J  Ihk 

Vfanpi  ",  fi:r»n  Mr.#-. 


8.  Right-click  on  Heat  Transfer  and  left-click  to  add  Heat  Source.  For  General 
source  Q,  Select  User  defined  and  Put  “f’,  which  is  the  heat  source  defined  from  the 
Variable. 

J  i  D^h;^  ■  h  Heat  Source 

--  Globa?  Definition  =■ 
j  —  Model  1  fmod-U 
^  Definitions 
j  -  Geometry! 


f\ Interval  I  (K) 

^  F  orm  Urmxi  (fin} 
k  !§;  Maten-ali 

:  #  Self  Defined 
s  Z  Hent  Transfer  (ht) 

•*  Hent  T ranker  an  Solids  1 

*  Th^rTTiisI  Insulfltior.  1 
l*  Mifll  VflbH  t 
Hrnt  Source! 
j  1^  M«h  1 

£i^ 

101 

i  (fg  Study  1 

-  Slc|j  1:  from  0  LO  ptiud  ifiC 
■  Solver  CohTitjurdliofi& 
i  Q  Results, 


»  Heat  Source 

o  General  source 

0 

User  defined 

T 

Wm 


Linear  iourt,  t 

Q-t^'T 


9.  In  Initials  Values,  make  Temperature  equal  0. 
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SEttings  .  1  Model  Library 


s  Initial  Values 


QomaEns 

Selection:  ^11  domains 


del  Guilder 

1-D.mph  (root} 

E  Global Definrtaans 
-  Model  1  (modi) 

E  Definitions 
#  GecmeifYl 
^  Interval  1  (tIJ 
Form  Union  ffinj 
d  Sj  Materials 

$  Self  Defined 
d  |  Heat  Transfer  (ht) 

Heat  Transfer  in  Solids  I 
S  Thermal  Insulation  1 
-  Initial  Values  L 
Heat  Sou  reel 
j  ®  Mesh  I 

jjs* 

l  Edge  1 
i^Ti  !>tudyl 

[u  > Step  1:  from  0  to  period  sec 

h  ?•*  SnEvpr  fnnfioiifaticn=:  ^  ^  ^ 

10.  Right-click  on  Mesh  and  select  Edge,  in  Edge  under  Element  Size,  select  Custom  and 

make  Maximum  element  size  equal  1/200.  Then  select  the  build  all  button  to  build 
the  mesh. 


T  Initial:  YaFues 
Temperature: 


1  0[KJ 

Surface  rad iosity: 


Model  Builder 


!  i  1-D.mph  (root) 

D-  E  Global  Definitions 
a  -  Model  1  [modi) 

>  E  Definitions 
a  ''  Geometry  1 

Interval  1  fil) 

Form  Union  fjfinj 
a  ®  Materials 

[>  E§3  Self  Defined 
a  §  Heat  Transfer  (ht) 

l<  HeatTransfer  in  Solids  1 
l<  Thermal  Insulation  1 
??  Initial  Values  1 
Heat  Source  1 
a  @  Mesh  1 

M  Size 
E3  Edsel 
al Sizel 

a  (f|  Study  1 

h\  >  Step  1:  from  0  to  period  sec 

r.  Pa  - ,4-1 


Settings  [00  Model  Library- 


a  i  I  Q 


^  Size 


Geometric  Scope 


G  eo  m  etri  c  entity  I  evel :  Enti  re  g  eo  m  etry 


EEement  Size 


■ . 1  Predefined  |  Extremely  frne 
o  Custom 


^  Element  Size  Parameters 
171  Maximum  element  size: 


1/200 

□  Maximum  element  growth  rate: 


1.1 


O  Resolution  of  narrow  regions: 

IT" 


11.  Under  Study  in  Step,  select  Range  button 


,  under  Entry  method,  select 


Number  of  values,  start  from  0  and  stop  at  “period”  which  is  defined  in  the  parameters. 
Put  Number  of  values  to  be  20. 
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del  Builder 
1-D.mph  (root) 

F  Global  Definitions 
—  Model  1  (modi) 

Definitions 
A  Geometry  1 

Interval  1  (il) 

Form  Union  (fin) 

5§3  Materials 

£§3  Self  Defined 
^  H eat  Tran sf  er  (h t) 

H  eat  Tra  n  sf  er  i  n  Sol  i  d  s  1 
^  Thermal  Insulation  1 
Sf  Initial  Values  1 
Heat  Source  1 
^  Mesh  1 
M  Size 

E3  Edgel 
M  Size! 
m  Study  1 

|yv  >  Step  1 :  fro  m  0  to  p  eri  o  d  sec 
fhv.  Solver  Configurations 
pT]  Solver  1 

Compile  Equations:  from  0  to 
u.v.w  Dependent  Variables  1 
|/V.  Time- Depen  dent  Solver  1 
El  Direct 
Jm,  Advanced 
■p  Fully  Coupled  1 
Di recti 

Q  Results 

||f|  Data  Sets 

Derived  Values 
gU  Tables 
h=t  ID  Plot  G  roup  1 
If]  Report 

@1  Player  1 


do  del  Library" 


?*£  Settings 


1a.  Time  Dependent 


a  =  n 


▼  Study  Settlings 


Ti  m  es:  rang  e(0f  (p  eri  o  d  [1/s]  -0)/19,  p  eri  o  d  [1/s]) 


s  l-J 


Range 


Entry  method: 
Start: 

Stop: 

Number  of  va  I  u  es : 


!  N  u  m  b  er  of  va  I  u  es  ^ 


period  [1/s] 


20 


Function  to  apply  to  all  values:  |  None 

j  Replace  |  Add  |  |  Cancel 


Mesh:  Mesh  1 


Physics  Selection 
Physics  interfaces: 


Heat  transfer  (ht) 


Use  in  this  study 


Discretization:  Physics  settings 


12.  Under  Time-Dependent  Solver,  make  Relative  tolerance  to  be  1.0e-6.  Left  Click 
on  to  compute. 


s>idi«l  Ri  lil.'l  hi 
■  l-[>.rpph  (met) 

—  Eilobat  DcflHrtiOns 
k+odeJi  1 

—  Prfipi(M?P5- 
e\  Ecornc-tjy  1 

Interval  1  (lJJ 
tffljl  Form  Union  (fm) 

Serf  Defined 
5  1  lent  Transfer 

I  Ht  ill  1  rrt-rit.1V  i  in  ‘tulidhli 
A  I  hernial  Insulation  1 

II  Initial  Values  1 

Ht-  rll  ^l?l  II  I  h-  1 

.5  Hcshl 

^  5i=* 

;  F I  I IJ  *•  1 

J 

xi  Study  1 

.■"^iHfi  I  :  frcnii  fl  I c :  | : r i : :  id  ‘.nr 

■" ..  Solver  Contigurabons 
Solver  1 

fi:mpilH  Fc|n>il,inn'^  frmmfl 

UcpL-i  d L n r  Vane  bios  1 

Time-Dec  endent  Solver  1 

El  Direct  ~ 

AibaMld 
»  fully  Coupled  1 
E3  Direct! 


•^H+ri*tg«, 


Mnrdri  I  iln^ry 


■r  -  I  U 


:l  limc-lJcpondont  Uolver 

T  General 

TirriOt: 


ra  itgcfd.  iptTi u  d|".  ft  |  J  i/IO.  pei  iudri/sU 

5  [uj 

1  IIh-B 

►  Absolute  Toleranze 

►  lime  Stepping 

k  Ki.-^ul  tv  WIiiIl-  S-ulwimj 
k  Oulpul 

►  Advanced 
b-  Lo  g 


13.  Finally,  under  Results,  select  Line  Graph  under  ID  Plot  Group  and  select  time  to  be 
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1  only,  we  should  be  able  to  see  the  result  like  Figure  9. 


14.  Further  data  analysis  can  be  done  under  Report,  such  as  generate  a  movie  and  export 
data  and  further  compare  results  with  MATLAB. 
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APPENDIX  D.  CRANK-NICOLSON  CODE  FOR  2-D  SIMULATION 


N=128;  %  number  of 
M=128;  %  number  of 
dx=l/N; 
dy = 1 / M ; 

x— ( [ 1 : N] '-0.5) *dx; 
y=( [ 1 : M] '-0.5) *dy; 
[xx, yy] =meshgrid (x, 
u=zeros (M*N, 1) ; 
d=  0 . 0  2  ; 
a=l ; 

u0=zeros (M*N, 1) ; 
u=u0 ; 


numerical  points  in  x  direction 
numerical  points  in  y  direction 
%  spatial  step 

%  numerical  grid 

y) ; 

%  solution  at  current  time 
%  radius  of  Gaussian  source 
%  magnitude  of  Gaussian  source 
%  initial  value 


dt=0 . 5*dx; 
m=  1 1 ; 

t= [ 0 : m] *dt ; 


%  time  step 
number  of  time  steps 


rl=dt/ dxA2 ; 
r2=dt/ dyA2 ; 


A  =  sparse ( [ ] , [ ] , [ ] , N*M, N*M) ; 
B  =  sparse ( [ ] , [ ] , [ ] , N*M, N*M) ; 
%%  forming  the  matrices 
%  the  bottom  boundary 
k=l ; 

A(k, k) =l+rl/2+r2/2; 

A ( kf k+1) =-rl/2; 

A ( k r  k+N) =-t2/2; 

B ( kf k) =l-rl/2-r2/2; 

B (k, k+1 ) =rl/2 ; 

B ( kf k+N) =r2 / 2 ; 

o 

o 


k=N; 

A(k,k)=l+rl/2+r2/2; 
A (k, k-1 ) =-rl/2 ; 

A (k, k+N) =-r2/2 ; 
B(k,k)=l-rl/2-r2/2; 
B (k, k-1) =rl/2; 

B (k, k+N) =r2/2 ; 

o 

o 


for  i=2 : (N-l ) , 

A(i,i)=l+rl+r2/2; 
A  (if  i-1 )  =*— rl/2 ; 

A ( i , i+1 ) =-rl/2; 

A ( i , i+N) =-r2 / 2 ; 

B ( i , i) =l-rl-r2 /2 ; 
B (if i-1 ) =rl/2 ; 

B (i r  i+1 ) =rl/2 ; 

B ( i f  i+N) =t2/2 ; 

end 

%  the  middle  layers 
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for  j=2 : (M-l) , 

for  i=2 : (N-l) , 

k=i+ ( j-1) *N; 

A (k, k) =l+rl+r2; 

A (k, k-1) =-rl/2; 

A (k, k+1) =-rl/2; 

A (k, k+N) =-r2/2 ; 

A (k, k-N) =-r2/2 ; 

B (k, k) =l-rl-r2; 

B (k, k-1) =rl/2; 

B (k, k+1) =rl/2; 

B (k,  k+N) =r2/2 ; 

B (k,  k-N) =r2/2 ; 

end 

end 

%  left  boundary 
for  j=2 : (M-l ) , 

k=l+ (j-1) *N; 

A (k, k) =l+rl/2+r2; 
A (k, k+1) =-rl/2; 

A (k, k+N) =-r2/2 ; 

A (k, k-N) =-r2/2 ; 

B (k, k) =l-rl/2-r2; 
B (k, k+1) =rl/2; 

B (k,  k+N) =r2/2 ; 

B (k,  k-N) =r2/2 ; 

end 

%  the  right  boundary 
for  j=2 : (M-l) , 

k=N+ (j-1) *N; 
A(k,k)=l+rl/2+r2; 
A (k, k-1) =-rl/2; 

A (k, k+N) =-r2/2 ; 

A (k, k-N) =-r2/2 ; 

B (k, k) =l-rl/2-r2; 
B (k, k-1) =rl/2; 

B (k,  k+N) =r2/2 ; 

B (k,  k-N) =r2/2 ; 

end 

%  the  upper  boundary 
k=l+N* (M-l) ; 
A(k,k)=l+rl/2+r2/2; 

A (k, k+1 ) =-rl/2 ; 

A (k, k-N) =-r2/2 ; 
B(k,k)=l-rl/2-r2/2; 

B (k, k+1) =rl/2; 

B (k, k-N) =r2/2 ; 

o 

o 

k=N+N* (M-l) ; 
A(k,k)=l+rl/2+r2/2; 

A (k, k-1 ) =-rl/2 ; 

A (k, k-N) =-r2/2 ; 

B (k, k) =l-rl/2-r2/2; 

B (k, k-1) =rl/2; 

B (k, k-N) =r2/2 ; 
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for  i=2 : (N-l) , 

k=i+N* (M-l ) ; 

A(k, k) =l+r l+r2 / 2 ; 

A(k,k-l)=-rl/2; 

A ( k, k+1) =-rl/2; 

A(k, k-N) =-r2 / 2 ; 

B ( k , k) =l-rl-r2/2; 

B(k,k-l)=rl/2; 

B ( k, k+1 ) =rl/2 ; 

B ( k, k-N) =r2/2; 

end 

%%  build  the  vector  f 
f=zeros (N*M,m) ; 
f l  =  zeros (N,  M)  ; 
for  k=l : m, 

xc=0 .5+0. 25* cos ( 10*2*pi*t (k)); 
yc=0. 5+0. 25*sin( 10*2 *pi*t (k)  )  ; 

f f=a*exp ( - (xx-xc) . A2/ (2*dA2) - (yy-yc) . A2/ (2*dA2) ) / (2*pi*dA2) 
f ( : r  k) =re shape ( f f , N*M, 1 ) ; 

end 

figure (2 ) 

drawnow 

for  k=l : (m-l ) , 

b=B*u+dt* (f ( : , k) +f  (  : , k+1) ) /2; 

u=A\b; 

pause  (0.2) 

temp=reshape (u,M,N) ; 
h=surf (xx, yy, temp) 

set (h,  1 edgecolor '  ,  ' none  1  ,  ' f acecolor '  ,  ' interp ' ) ; 

%view ( 2 ) 
drawnow 

end 
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APPENDIX  E.  FFT  CODE  FOR  2-D  SIMULATION 


clear 

o 

o 


N=12  8 ; 

dx=l/N; 

x= [ 0 : N ]  1  * dx ; 


number  of  numerical  points  in  [0,  1] 

h  spatial  step 
b  numerical  grid 


[xa,  ya] =meshgrid (x, x) ; 


u=zeros (N+l, N+l) ; 
d=  0 . 0  2  ; 
a=l ; 

u0=zeros (N+l, N+l) ; 
u=u0 ; 


%  solution  at  current  time,  u(j)=u(x(j)) 
radius  of  Gaussian  source 
magnitude  of  Gaussian  source 
%  initial  value 


dt=dx/ 8 ; 
T=12/128 ; 
m=T/dt ; 
t= [ 0 : m] *dt ; 


h  time  step 
number  of  time  steps 


o 

o 

%  going  to  the  coefficients  of  cosin  expansion 
w= [ u ;  u(N:— 1:2,  :) ] ; 
w= [w,w(:,N:-l:2) ] ; 
z=fft2 (w) ; 

cu=real (z (1 : N+l, 1 : N+l ) ) /NA2  ; 

Na= [0 :N] ; 

[Nx,  Ny ] =meshgrid (Na,  Na) / 
r=(Nx*pi) .A2+(Ny^pi) .A2; 
r ( 1 , 1 ) =1 ; 

hl= (1-exp (-r*dt) ) . /r; 
hi (1,1) =dt ; 

h2=(l- (r*dt+l) . *exp ( -r*dt ) ) ./r.A2; 
h2 ( 1 , 1 ) =0 . 5*dt A2 / 
r (1, 1) =0; 

o 

o 


x0=0.5  +  0.25^cos ( 10^2^pi^0 )  ; 

y0=0 . 5+0 . 25^sin ( 10*2*pi*0 ) ;  %  location  of  Gaussian  source  at  tA{k-l} 

f0=a*exp (- ( (xa-xO) . A2+ (ya-yO) . A2) / (2*dA2) ) / (2^pi^dA2 ) ;  %  Gaussian 

source  at  tA{k-l} 
w= [ f  0 ;  f 0 (N : -1 : 2 ,  : ) ] ; 
w  =  [  w ,  w  (  :  ,  N  :  - 1 :  2 )  ]  ; 
z=fft2 (w) ; 

cf 0=real (z ( 1 : N+l , 1 : N+l ) ) /NA2; 

o 

o 


surf (xa,  ya,  u,  ' edgecolor ’  ,  ' none '  ,  ' f acecolor ’  ,  ' interp ' ) 

axis ( [ 0 ,  1,  0,  1,  -0.04,0.4]) 

caxis  (  [ 0 ,  0.4]) 

view ( 3 ) 

drawnow 

for  k=l : m, 

xl=0.5+0.25*cos ( 10*2*pi*t (k+1) ) ; 

yl=0 . 5+0 . 25*sin (10^2^pi^t (k+1) ) ;  %  location  of  Gaussian  source  at  tAk 

f l=a*exp (- ( (xa-xl) . A2+ (ya-yl) . A2) / (2*dA2) ) / (2^pi^dA2) ;  %  Gaussian 

source  at  tAk 
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w= [ f 1 ;  fl (N:-l:2, : ) ] ; 
w = [w,w(:,N:-l:2) ] ; 
z=fft2 (w) ; 

cf l=real (z ( 1 : N+l , 1 : N+l ) ) /NA2; 

%  update  cu 

cu=cu  .  *exp  (-r*dt)  +cfl.*hl+  (cfO-cfl)  /dt . *h2; 
cf 0=cf 1 ; 
k 

%  going  back  to  the  function 
z=NA2  * [cu;  cu (N:-l :2,  : ) ] ; 

Z= [z, z ( : ,N: -1 : 2) ] ; 
w=ifft2 (z) ; 

unreal (w(l:N+l,l:N+l) ) ; 

%pause (0.2) 

surf (xa,  ya,  u,  ' edgecolor ' ,  ’ none 1 ,  1 f acecolor 1 , ’ interp ' ) 
axis  (  [  0 ,  1,  0,  1,  -0.04,0.4] ) 
caxis  (  [ 0 ,  0.4]) 
view ( 2 ) 
drawnow 
end 
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APPENDIX  F.  COMSOL  CODE  FOR  2-D  SIMULATION 


2.  Open  COMSOL  4.0a  with  2D  and  hit 


ft*  Model  Wizard  H  Model  Library  Material  Browser 

=  □ 

Select  Space  Dimension 

p  ^  ^ 

C  3D 

2D  axisymmetric 

o  2D 

ID  axisy  m  metric 

ID 

0D 

2.  In  Heat  Transfer,  select  Heat  Transfer  in  Solids  (ht)  and  right-click  Add  selected  then 
rit  . 

Add  Physios 

0  Recently  Used 

t>  AC/DC 
|>  i=3  >■  Acoustics 

!.-  Chemical  Species  Transport 

O  ij _ ff  Electrochemistry 

l>  Fluid  Flow 
^  \S\  Heat  Transfer 

^  Heat  Transfer  in  Solids  (ht) 

SEjl  Heat  Transfer  in  Fluids  (f  ^  Add  Selected 

c3jjl  Heat  Transfer  in  Porous  rvrtitarencm^'' - 

s  Heat  Transfer  with  Surface-to-Surface  R.adiation  (ht] 

3  i  o  h  eat  Tra  n  sf  er  (ht] 

#  Surface-to-Surface  R.adiation  (rad] 

Joule  Heating  (jh] 
f:  Q  P I  a  sm  a 
:  Au  Mathematics 


^  X 

Selected  physics 

^  Heat  Transfer  (ht] 

3.  In  Preset  Studies,  select  Time  Dependent  and  hit  Finish  ^  . 

Model  Wizard  H®  Model  Library:  ^  E3 

Select  Study  Type  <=>  i=>  [gj 


Studies 


Selected  physics 

±  H  eat  Tra  n  sf  er  (ht) 
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4.  Under  Model  Builder,  right-click  on  Global  Definitions  and  left-click  to  select 
Parameters,  input  Parameters  as  following: 

^|7  Mcdel  Bl  der  '  Setting:  [jjl]  Model  Libr^n1 :  MEterial  Grcwser 


■  i  2-D.nph  raj 
E  Globe  Definitions 
Parameters 
Mocel  Wl! 
£  Study  1 
Q  huis 


p.  Parameters 


Parameters 


NamE 

B:pnessic.n 

Value 

Description 

d 

0jfl2[m] 

0jO2m 

rac:  js  of  Gaussian  Souris 

period 

iw 

Is 

rotating  peried  of  neat  source 

3 

025[m] 

0.25  m 

rac  js  of  rc-atirg  tractor/  in  k  direction 

b 

Q.25[m] 

0.25  m 

radius  of  rotating  tr-  ectory  in  y  direct  on 

ifl 

05[m] 

05  m 

center  of  rotert  on  ir  :<  direction 

V3 

Q5[m] 

15  m 

center  of  lutat  on  ir  y  direction 

n 

IXJEfl  W/m2 

intensity  of  lie.:",  source 

5.  Right-click  on  Definitions  and  left-click  to  select  Variables  input  Variables  as 
following: 


Settings 


Mcdf!  Library  Material  Browser 


Entire  mode 


1;  Model  Builder  ^  D  B 
I  2-D-mpfi  fnppfj 
21  Global  Definitions 
Model  I  (modi) 

S  Definitions 
3=  Variables  I 

Boundary  Systcr 
^4  ViwI 
Geometry'  1 
i@r  Materials 
1  HeatTran-sfer  |7if.) 

@  Mesh  1 
Jiii  5tudy  L 
G  Result? 

6.  Right-click  on  Geometry  and  left-click  square  and  make  side  length  1  m.  Click 


*=  Variates 

Geometric  Scop*: 
Geometric  entity  level: 

▼  Variables 


Name 

Expression 

Unit 

KC 

sO  +■  e-’co  s[2*p  Wperiod 

m 

yc 

/)+  h J  sin(2*pi"l/p  eriodj 

m 

q 

[0/(2*pi*d  Al  [l/mA2])*&p£-ttJ(-Kcl  A2-  <y-yc  JA2)/(2*d  AZ)  J 

Wm* 

build  all  button  to  create  a  square. 


Mnrli-’  Risilrdhi  tt-  d  ^ 

+.1  2-[>  rnpli  I’rcritJ 

=  l-ilufcjd  Definitions, 
Mod*.  1  i'bfOo4.i 
ez,  Defipi+iont 

C-i  iji  131  I  iy  1 

_  |  SqLarcl  zql.f 
Form  Union  iYTn 
S'"  Mill  h  M.ll  . 

4  Heat  I  ransfer  (hi) 
£5)  M«h  t 
p---^  ^inilyl 
E3  Results 


'n fc-- 1- i r i ij >•  I]  Mrrrtnl  |  Mh^mhI  Rmk-jKM 

CJ  1 

Ir  1  U  “  n 

□  Square 

t-  Object  Type 
Type;  |  Solid 
■—  Size 

Side  length:  1 


J^;c.  Cc- frier 


k:  a 

V" 

-v  Rotation  Ang  to 


KotjfciOri:  0 

*  Lavzfi 


deg 
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7.  Right-click  on  Materials  and  left-click  on  add  Materials,  add  Steel  AISI  4340. 
Right-click  on  the  “square”  in  Graphics  and  left-click  on  it  to  select.  Make  sure  that  “1” 
is  under  the  selection. 


Model  Builder 

d  *  2-D,mph 

=  Global  Definitions 
4  |  '  Model  l  (modi) 

\>  S  Definitions 
t?  /\  Geometry  1 
4  y  Msterialr 


#  Sled 

0  m  Self 

|  HeatTn 
0  &  Mesh  1 
■  ^  Study  1 


Settings 


Model  Lit 


§  Materials 


Material 

Open  Material  Browser 
Dynamic  Help  FI 


©  Results 


■■■?>.  hj  J.' 

■  n  Z-a-rn  r.iT. 

T  n  LTrd  :tj 

1  "Vi."  Lfeiri 

h  ferfU  .  r. 

i 

3  IM-irirH 

rr  Mm3 


3  |  *  1  'id>v  Mwrll  hn-?  #  Wjlr>  A  few 
«  HrftttJ 
CiwplMlIiipi; 

■J.TTd*-;  Ti+I  Kvd  [Uh  k  ■ 

;:■*  J?r  N:«js 


\  + 

r  - 
* 


>V--kT :  t+rvnY 
.IriTTijr*  r  -.'rid  i 
he  *  MezI-stu 

■h  -|i-h>  hi 

3hMvU 


I  NlW'. 


.'r:dt.  |  |Uf* 

-h  -V-i  >,  *  .  T|i-.ir  l'|i  J.'J  |Lli5'Li 

"ifir.  r>H  TK^^Tip]1  nj-^'1 

■4tx  kt : >:t.+.  I  ^ Jh'.' A-- It "*■ 


L  Li  -  ■  ’  ■  \  ■  *  ~  = 


M-  EZ  i.i.j  14  MJi  Wt  I.'  Wt  J.J 


Ilt  i  :-|  =  '  :  1r  j-Cij 

IWiV.ll’JUl: 


8.  Right-click  on  Heat  Transfer  and  left-click  to  add  Heat  Source.  For  General 
source  Q,  Select  User  defined  and  Put  “q”,  which  is  the  heat  source  defined  from  the 
Variable. 
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i.  MudtJ  Bui Wfcr 
j  lJP  2-Dumpti  frcdtJ 

=  Globs  i  Otfi  mtiori  s 
m  Model  L  iVroa'i; 

Definrti-cins 
'/'•  Grom-rdry  1 
I  ®  Material-? 
j  ^  Hfcal  Tr-dsiiPer  ffrU 

*  H  cot  T ranker  Ln 
"  Thermal  2njulati 

*  Initial  Vdlutt  1 
Heat  Source  1 

Me?h  1 
STEJ  ^  U  1 
r-  Si.1  Results 


StLlinyt  O  Model  Library  Material  bru-LVLti 

Heal  ^ourw 

EhwqalpK 


a  “  ^ 


Selection  Manual 


w  Meat  Source 

O  Otr' fdl  vL'Uir.  t 


I?  User  defined 


W/mJ 


linear  source 

Q=rfcT 

9.  In  Initials  Values,  make  Temperature  equal  0.  We  assume  there  is  no  temperature  rise 
in  the  beginning. 


".  Mode' Builder 
j  i  2-D.mph  (root) 

.7  Gloh.il  rVrir.ihons 
4  \  Model  1 1'modJJ 
^  DEfinrtions 
/ \  Geometry! 
Malt-lidk 

#  ^  H  cot  T ran  sf  er  ft  t) 

LjJ  H  Eat  T ran  &f  er  in 
^  ThErmal  ]ns-ulab 
^  Initial 'idluei- 1 
I2  Nest  Source  1 
Q  Mesh  I 
(§ti  Study  t 
Q  Retultt 


Settings  ill]  Model  Library 

Initial  Values 


Material  Broker 


a 


Domain* 

Selection: 


\  4= 

ra  = 

X 


»  Initial  Value? 
Temperature 

t  0[iq 

Sumaie  -adiccib/: 


J  0 


0.  Right-click  on  Mesh  and  select  Mapped.  Right-click  on  Mapped  and  select  Edge 
Groups.  Please  select  domain  and  add  each  edge  group.  Under  size,  select  Custom  and 
make  Maximum  element  size  equal  1/64  or  any  1/2AN  where  N  is  an  integer  to  compare 


with  FFT  2  method  in  MATFAB.  Then  select  build  all  button 


to  build  the  mesh. 
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Put  Number  of  values  to  be  3 1 . 
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pfl Model  Builder 
■  i  2-D.rnph  (root) 

B  Globa!  Definitions 
Pi  Para  meters 
l!  Model  1  (modi) 

S  Definitions 
A  Geometry  1 
ij$  Materials 
|  Heat  Transfer  fjhy 
^  Heat  Transfer  in  Solids  1 
co  Thermal  Insulation  1 
L-'  Initial  Values  1 
Heat  Source  1 
H  Mesh  1 
^  Study  1 

1/\  >  Step  1:  Time  Dependent 
[(Yp  Solver  Configurations 
Q  Results 


Settings  [HO  Model  Library  Si*  Material  Browser  dL  a  B 

Time  Dependent 

▼  Study  Settings 

Times;  rang  etO,  (peri  o  d  [1/s]  ^3)/30,period[l/sJ) 


y 


Range 


Entry  method; 

Start: 

Stop: 

Number  of  values:  31 

Function  to  apply  to  all  values: 


Number  of  values  ^ 


periodil/s] 


None 


Replace 


Add 


Cancel 


2.  Under  Time-Dependent  Solver,  make  Relative  tolerance  to  be  1.0e-5.  Left-click 


on  to  compute. 


1  [I  I  Model  Builder 

■  »  2-D.mph  i>OL'L.j 

Global  Definitions 
Pi  Parameters- 
I  Modol  1 

=  Definitions 
lj  eo  m  etry  1 
gfcK  MdLendlir 
^  Hrnt  T  rnn^for  f’htj 

Heat  Transfer  in  Solid s  1 
,Ji  Thermal  Insulation  1 
^  LrnLidl  Values-1 
Host  '>l-j lj rr c  1 
Mesh  1 
d£iil  Stu  dy  1 

|  '! ,  >  Step  1:  T  itti e  Depend enL 
Tr n,  ■^nlvr-r  C.  n-nfiqLjrntmn-q: 

H  Solver  1 

jt-L1!  Compile  Equations:  Tin 
D ep en dtnL  Variables  1 
|.’^:>  Timr  D op on dent  ^,nlvo 
|^|  Direct 
2u  Ad va  n  c  ed 
■r  '■  F-ully  Coupled  1 
E5I  Dirortl 

L.  i  Results 


Set  Lings  II I U  Model  Library  35"  Mete  riel  BrowLtn 


Rv  lime  Dependent 

Times:  renget,Ortpiericid[l/E,]-03/'JOi,peTi'odIl  s 


Relative  Lo  I  era  nee:  l.QeO 


►  Absolute  Tolerance 

►  Tim  e  St  ep  pi  ng 

►  Resul  ts  WhiTe  Solving 

►  Output 

►  Advanced 

►  L«j 


4.  Finally,  under  Results,  select  Surface  under  2D  Plot  Group  and  select  time  to  be  1 
only,  we  should  be  able  to  see  the  result  like  Figure  25. 
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,  .uai 


File  Edit  Options  Help 

PBBfl □ a i a  ' 


t7  Model  Builder  7°f 

V»  2-D.mph  (root) 

£  Global  Definitions 
Pi  Parameters 
ft  Model  1  (modi) 

£  Definitions 
A  Geometry  1 
$  Materials 
J  Heat  Transfer  (by 
f  Heat  Transfer  in  Solids  1 
5  Thermal  Insulation  1 
Initial  Values  1 
*  Heat  Source  1 
^  Mesh  1 
m  Study  1 

]£,  >  Step  1:  Time  Dependent 
[^,  Solver  Configurations 
§]  Solver  1 

fjj*  Compile  Equations:  Tin 
u,v,w  Dependent  Variables  1 
^  Time-Dependent  Solve 
[S  Direct 
Advanced 
Fully  Coupled  1 
S  Direct  1 

Q  Results 
i  DataSets 
Views 

[«  Derived  Values 
H  Tables 
□  2D  Plot  Group  1 
Q  Surface  1 
D=1  Report 


■  Q  fc  Q  g|  ^  H  B _ _ 

Settings  Hi  Model  Library  $  Material  Browser  a  3  ,;l  G 


'I  a 


Q  Surface 

▼  Data 


Data  set  |  Solution  1 

Time  [l 


3B 


▼  Expression 
Expression: 


E 


□  Description: 


[Temperature 


►  Range 

▼  Coloring  and  Style 
Coloring:  [Co 

Color  table:  [to 

12)  Color  legend 
□  Wireframe 

►  Quality 

►  Inherit  Style 


$  Q.  &  EE  1 1®  “  Q 


Surface:  Temperature  (K) 


14.  Further  data  analysis  can  be  done  under  Report,  such  as  generate  a  movie  and  export 
data  and  further  compare  results  with  MATLAB. 
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APPENDIX  G.  COMSOL  CODE  FOR  3-D  SIMULATION 


1 .  Open  COMSOL  4.0a  with  3D  and  hit 

*  ■*  Mcdd  Wizard  —  -■ 


Select  Space  Dimension 

e  3D 

2D  axisymnietfit 

©2D 

ID  axisynursttric 

©ID 

DO 

2.  In  Heat  Transfer,  select  Heat  Transfer  in  Solids  (ht)  and  right-click  Add  selected  then 
rit  . 

Add  Physics 

0  Recently  Used 

t>  AC/DC 
|>  <]->■  Acoustics 

!,  Chemical  Species  Transport 

O  ij _ ff  Electrochemistry 

l>  Fluid  Flow 
^  \S\  Heat  Transfer 

^  Heat  Transfer  in  Solids  (ht) 

SEjl  Heat  Transfer  in  Fluids  (f  ^  Add  Selected 

c3jjl  Heat  Transfer  in  Porous  rvrtitarencm^'' - 

s  Heat  Transfer  with  Surface-to-Surface  R.adiation  (ht] 

3  i  o  h  eat  Tra  n  sf  er  (ht] 

#  Surface-to-Surface  R.adiation  (rad] 

Joule  Heating  (jh] 
f:  Q  P I  a  sm  a 
:  Au  Mathematics 


*  x 


Selected  physics 

^  Heat  Transfer  (ht] 

3.  In  preset  Studies,  select  Time  Dependent  and  hit  Finish  ^  . 

J  Model  Wizard  |®]  Model  Library;  ^  E3 

Select  Study  Type  <=>  ■=>  |^~ | 


Studies 


Selected  physics 

±  H  eat  Tra  n  sf  er  (ht) 
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4.  Under  Model  Builder,  right-click  on  Global  Definitions  and  left-click  to  select 
Parameters,  input  Parameters  as  following: 


'  Mode!  Builder 

j  ■  >  3-B_mesti.mph  iroct} 
a  F  Global Def in rtiens 
pi  Parameter 
i  Moddl  fm«iliJi 
Study  1 
O  Results 


Settings 

r,  Parameters 

t  parameter? 


Name 

Eipre^ion 

VftlUE 

|  Description 

d 

0.02  [m] 

P.G2  m 

oenod 

tW 

Is 

a 

O.Z5|m] 

D.Z5  m 

radius  of  f-otitin...y  in x. direction 

b 

0t2S[m] 

0.25  m 

radius  of  rptatin..,y  in  y  direction 

xD 

0.5[m| 

D.5m 

CEntEr  of  rotation  in  a  d  rection 

0.5[m] 

D^m 

center  of  rotation  in  y  direction 

a 

5.0eS[W/m  A2] 

500000,0  W/ma 

5.  Right-click  on  Definitions  and  left-click  to  select  Variables.  Input  Variables  as 
following: 


Entire  medei 


|  ,  Model  Builder 

■  i  3-&_mBh,mph|VoDt 
=  Global  Definitions 
i  Mod  e!  1  \/7Jodl) 

=  Definitions 
a=  Variable^ 

J'  Boundary 
^  Viei^'l 
Geometry! 

$  Materials 
|  Hfcst  Transfer  (fttj 
©  Mesh  L 
5^  Study  1 
G  fiespjlt& 

6.  Right-click  on  Geometry  and  left-click  “Work  Plane”  and  input  1  in  z  coordinate. 


Settings 

Variables 

Geometric  Scope 
Geometric  entity  level: 

t  Variables 


Nsme 

impression 

|  Unit  Des.Jon  | 

3& 

KO+a^ccsI^piVperind] 

m 

y* 

yO+b^sinCl'pl^ptricHf) 

m 

q 

U/(2‘p  i1  d  A2:  [i/m  A2])  fc&:p(-((x-jcc)  *2-  (y-yc)  A2]/(2“d  *2)) 

Mm* 

Model  Builder 

8  1  Untitlp-d.mph  (nxtj 

r=  Global  Definitions- 
a  s  Model 

==  Definitions 
a  _  Geometry] 

9^  ■>  Form 
m  Materials 
^  Heat  Trans 
&  Mesh  1 
P  <*C=j  Study  1 

|>  i'i  RtbLlltb. 


Build  All 

Import 

Bio^k 
Curie 
Cylinder 
Sphere 

More  Primitives 
Work  Plane 


Settings 

Geometr 

Geometry  Si 

I  Ini4-r 

F& 


m - 

C 

r- 

t 


90 


7.  Right-click  on  Geometry  under  Wok  Plane,  select  Circle.  We  need  to  build  up  two 
circles.  One  has  radius  0.32m  and  the  other  one  has  radius  0.18m,  both  are  centered  at 
x=0.5m  and  y=0.5m.  Right-click  on  Geometry  and  add  Extrude,  select  Reverse  direction 
and  input  distance  to  be  0.05m.  Right-click  on  Geometry  and  add  Block  with  Width, 

Depth  and  Height  all  equal  1 . 

4  Geometry  1 

4  .  ''j  >  Work  Plane  1  (vrfil} 


>  ^  View  2 
>  Form  Un 
affc  Materials 


[> 


•  "-jr  Mesh  1 

^  CE.  .^4,1 

First  circle: 


Hi 

Build  All 

PS 

+G? 

Import 

( *D 

Cirde 

*(=> 

Ellipse 

■*  Untitled. mph  {rixt} 

GluL-dl  DcPimLiuritr 
J  Model  1  fmodV 
=  Definitions 

-  A  Gr-orrtCtry  [ 

m  ..7^  >  Warfc  Plane  I  fwpl J 
j  Geometry 

Q  Ciftlel  frJj 
EQ  View  2 
•=■  Form  Union  ffinj 
Materials 

^  I  Seat  Transfer  fhy 
Mesh,  1 
risj^l  Study  1 
Results, 

Second  circle: 


i  Untitled. mph  (rootj 

dubel  Definitions 
J  V  Mo  del  1  Cmodi? 

Definitions 
-  A  Geometry  1 

j  r'  >  Work  PEane  1  {wpl) 
j,  Geometry 

O  Circlet  fcl) 
Circle  ?  (c2) 

I-  ri-H  View  2 
^  >  F-o-rm  Union  (fin) 
MdLt-ldli. 
i  HedL  1  ransrer 

^  Mesh  1 
I..  Study  ! 

Extrude: 


Circle 

▼  Object  Type 
Typ-r:  i^rilifl 
-w  Size 
finding  G.lfi 
Position 
Jin  sr:  j  f  rntrr 
W  0r5 

y.  0.5 
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\  Model  I  fi-ncwd-D 
“■  Definitions 
yk  Geometry  1 

-§r  Work  Plane  1  fyvpl.} 
Extrude  1  f'extJ.J 
=>■  Porm  Union  fftnj 
%£=:  Mjlcrk3h 
^  H  v^il  F  E^cc  L^KiHr  |>liy 
NTiri-r-hi  I 


5tudyl 
fiSl  Results 


■—■  General 
Work  plane:  j  wpl 
Input  objects: 

Kl  i  i  III  Print  iheIifkej  Stn  Fir 

I _ |  Keep  input  objects 

|Vl  Keep  cross -section  a  I  faces 

Distances  from  Work  Plane 

Distances  ( mi 
Q.n^i 


I V|  Reverse  direction 


Settings 

*t  u  u  ( o  ck 


■»■  Oliiftff  1  Ty|Mr 

Type:  Solid 

j-jiie!  SI  upL- 
widthi  1 

UElprh.li:  1 

M  iri  gKL  1 


►  Puviliun 


t  Axis 

A»ik  |  t"  Hrl**>iiriri 

V 

y*  o 

i 

I11  Ru  !<■  4  iui  ■  /Vi  il|I  w 
>  Layers 

he  geometry  should  be  built  like  this: 

b  □!  a  ^  -  I  si  <=i  sn  i«n  I  .j-  ■>  ^  tii  L:  I  r  <a  a  t -_l||  ■ 


■>  C 


Add  a  Block: 

T  ^  Model  Builder 


!  tippcndhf  qmph 

“  Global  Uetimtiona 
n.  K/lnrl**l  1  fmnrflj 
—  Definitions 

A  G"“"-1'y 1 

Work  P I  e  ne  1  fTvfj-tJ 

iZ--.,  Ejct  r  u  eJo  L  fpjrfjji 

'  Slock  1  ftW-X) 

>  3-orm  Union  fftsij 
??jp  M  Ml-  a>-ri  mI  k, 

-=^  Neat  Transfer  fihfJ 

l£ri5j  Mi^h  I 
■i'a.’  Stu  dry  ± 
m  Faults 
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7.  Right-click  on  Materials  and  left-click  on  add  Materials.  Add  Steel  AISI  4340. 
Right-click  on  the  geometry  in  Graphics  and  left-click  on  it  to  select.  Make  sure  that  all 
the  portions  are  selected.  “  1”  ,  “2”  and  “3”  should  be  under  the  selection. 


McirJpi  FfuilHl#r 
i  ‘1  Y-H  nFph.mp  n  fmrVy 
—  G  Libdl  Dt-H-iil  ii;  -i>- 
d  L  M  J  JlI  1  y.<  hLHj'J  1 

ir-  Cl ■_ T 1 1 1 1 C i *_  -ii 
■"  Ueofnrtiyl 
#  £t  Ma+tnilr 

.  t±i  Steel  r.N  t>IO 

■J  H  fji+  7  r.%  i  ■/*  r  (h  i'J 
MlHHh  ] 

-  «T>  SLudy  1 
d  i„jj  Rl-jullii 

■  UlfftS  itts 

LUrr-  ed  vnluts 
g*  T-sril« 

F-  Pint  firnt  ip  I 

IFj  Pint  Lirnup  J 

Prubv  1C-  FIliL  Giuu|i  3 
i  L-  I  F.'_pj'.''L 

Player  1 

m  Aniir-s;ion  1 


Mjvreriftl 


hi-iTiHl  ii  vnlity  Ihvi-I:  Dninnir 

beiertor-  A|1  dnnr-s-ns 


►  Material  Frogjcnfw 
w  Material  Loulmti 


|  PTDfKfty 

|  Njnr-e 

|  Value 

|  Unit 

1  l«(t  CSP-SC'ty  pnecsL^e 

tp 

f^L  ’KJJ 

J.-'lkq'K; 

C'enpf,- 

rho 

:m...  -Jj 

knf-^3 

THlfn-n^i  condiichwity 

k 

8.  Right-click  on  Heat  Transfer  and  left-click  to  add  Boundary  Heat  Source.  For 
Boundary  heat  source  Qb,  enter  “q”,  which  is  the  heat  source  defined  from  the  Variable. 
Right  click  on  selection  9  on  Graphics  and  left  click  on  it  to  select. 


*  3-B_mesh,mph  (root) 
=  Global  Definitions 


Boundary  Heat 


d  \  Model  1  (modi) 

=  Definitions 
/\  Geometry  I 
n  «  Materials 
d  ^  Heat  Transfer  (ht) 

Heat  Transfer  in  Solids  1 
%  Thermal  Insulation  1 
=■  Initial  Values  1 

Boundary  Heat  Source  1 
0  &  Mesh  1 
Study  1 
a  Q  Results 

Data  Sets 
Derived  Values 
Tables 


Boundaries 


Selection: 


Manual 


▼  Boundary  Heat  Source 


Boundary  heat  source: 
Ob  q 
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Flounclnr/  Heat  Spurr-e 
hunttaricf 

SduLu-:  Mjii-j 


&  1 


■  !  unphc? 


I  U  3  :*  - 1  ^  ^  3,  ■  &  S I 


I  •■  ••:•:  I  a  a  c[g|g  ^ 
IT-.; 


=i 

\  'fr 

t  - 

X 

ij 

1 

r--' 

\ 

|g 

V 

mi 

—  Boundary  1 1  eat  Split* 

R-»  ndnn  vn  '.rail- it 

ajj 

| 

ss 

Or,  q 

i 

j 

i '  ■ 


U*k 


t.i ' 


>LJ 

'  "&5 


n  c 


9.  In  Initials  Values,  make  Temperature  equal  0.  We  assume  there  is  no  temperature  rise 
in  the  beginning. 


*  3  H  eat  Transfer  jVifl 

L-  H^i  Tranter  in  SeNds! 
^  Thermal  IrroulaticMil 
Jnnal  Values- 1 
Gouns-ary  Heat  Source  1 
I-  ©  Mesh  1 
£■  1 
j  □  Results 
[■■  3:j  Data  Set; 

:  ^  Derived  Values 
:  !=  tables 
iTl  JO  ?  lot  Group  1 
!■■  (5  30  Plot  Group  2 
■  Probe  ID  Plot  Group  3 


E  “ 

4'  X 


T  Initial  Values 
Temperature 


T  D[K] 

SuFoce  radicraty: 


J  0 


WVm 


10.  Right-click  on  Mesh  and  select  Free  Tetrahedral,  repeat  this  process  three  times. 

>  6§J  Mesh-i- - - - U - 


Study  1 
*  Cl  Results 
:■']  Data 


Build  All 


hS  Oeriu 


.  J  Free  Tetrahedral 

*  lls  .  Swept 

Boundary  Layers 
l£j3DPI  More  Operations 

First  Free  Tetrahedral  mesh:  selection  2  is  picked  and  makes  the  predefined  Element  Size 
to  be  Extremely  Fine: 


Table 
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1~.  HjJ-J  ?j  ll 

J  ■  ■  1-P  ITF^l  n't!  VwjJ 

=  GcfctfBf'.iltFt 

*  k  h':: J  L  itj-j.1 

I  llriTjWK 

$  h'niih 
|  Het  Indt-  rr) 

d  l|l  K'rCll 

d  T-:.  Ik  jij^J 

J  FTJ' 

K  liwl.iyw 

Rul 

>  *  hrt,-- 

d  3  R«lI; 

.  I  Did  ld..i 
J  II**  M  k k> > 

■  TtiL. 

E  IDfl-  filial 

C  IPFrLuwZ 

■  ftrtcJDFu  SlitS 
d  Pl^::t 

_  Ffc-at 

n  a™  nii  it  . 


-Jiiza 

Gumtkkvi 

f-KWTFTV  FTTH  IdVd  I  riTTdl  T 
H.  J.I  Lr 


.  !  ’■  ~  ^  I  '  G  fcl-  L. 

~|  a  fl  o  an  -  ai  £l  r^.  ft  E3  ■£■  -  i-  '  ■  |  n  F-  E 


- 


'S  + 

tj  - 

-  K 


DjiulEin 

i  LVJH  rfl  {Ldwruk  fri 


r.d!l>T< 

►  I  If  TaH-d1  "dT.® .  ‘.in 


Si;-.. 

W0:\ 

■tr. *.*fx  /  isaO\  1 

l  X  i  ./ 

wfesZmFit  i 

N 


£ 

u- 


a: 


Second  Free  Tetrahedral  mesh:  selection  3  is  picked  and  makes  the  predefined  Element 
Size  to  be  Extremely  Coarse: 


i  "¥*l  I  i  M. 

m  i-  HidJI  rwM.i 

I  Nrirhm 

y, 

it  M^l  i  Jj 

\  I  Fll  llWVC *■  T 

f  i  .■  L-'ri’-i  I 

J^ 

J  lrtfr'-IV*dIl 

d  -■  friti  Gi.h 
.  j  Umil 
a  .  fo-i 

JCWK 
i  r  ^  - 
a  rj  td-.i- 
;J  [kuks. 

y  L4dV  El  VdlMd 
p  Td  I  r  : 

G.  3IT1t.ji:je  L 
l*i  J."  H-  IrajE  J 
fifh-  |n  Fk-d  r.v.^.~ 

ji  FI  hwL 

L  Fti.H-L 
■  |Vi  .-.-d...  i 


Litcnririr  'iE-jif 


■j K-rmv  itt.  krit  J e n i n 
Mexi  rAii.i 


■  f'lr.Wirl  :  Jj  LIT  I. EE 

Ljtvxti 


mmmwr 

/■i/,'<ruc\ . s 


/■ 


3 


v- -  -  r^.;  X’  "" 


Third  Free  Tetrahedral  mesh:  selection  1  is  picked  and  makes  the  predefined  Element 
Size  to  be  Extremely  Coarse. 


.'.Vifc  lu:ir 

*  i  r.i-^M^Hi  rM' 

'  ilaW  EHrfMI 
d  !  4s*fI1  rvsZ- 
“  Dd-tui, 

'A  fahviHr,  i 
t  Hilwk 

|  HelI  TijdJTdtH1 

d  1  Hi-J.  I 

J1™ 

j  TuL'jklpJ 

a  ■  LdVi  I  Yd 
,|  r.EITU' 

d  F.'.i . 

v.l  rwF 

UuhL 

a  iy  Fiij;, 

■ !  IUh'n 
I  'i  Lki  dir'.j  lii 

f3 

r,  1-Fix '-L-,  ^1 

Fll h£ JD  1  Ll  $ILL|.  3 

d  !  Fdaxl 

U  FL.-u  L 

B  i. F’d.1  L>  L 


I  I,  Q-  "1 


LietkIt:  drl  K'Il  tl  .enri 

t+flld-d 


:h  v*- 


L1C 


.•  In  ^  tS 

^■1'  .  i 

\i. 


d*  a  d-:  i+rv-:  Eiimd-Y 

arfin 

■  FrFicil  Sue  Fa-meIeif 


/-  ■  v-  \/  \  i 

V  #  Jp  k 


;--k 

/I 

^  i 


I 


u 


■in\K 


"■W 

d  : 


Then  select  build  all  button 


to  build  the  mesh. 
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,  under  Entry  method,  select 


1 1 .  Under  Study  in  Step,  select  Range  button  ' 

Number  of  values,  start  from  0  and  stop  at  1 .  Put  Number  of  values  to  be  5 1 . 
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12. 


Under  Time-Dependent  Solver,  make  Relative  tolerance  to  be  1.0e-5.  Left-click 

on  to  compute.  This  3D  problem  may  take  couple  hours  to  solve,  in  order  to  see  a 
quick  solution,  we  can  make  all  meshes  to  be  extremely  coarse. 
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2.  Finally,  under  Results,  select  Surface  under  3D  Plot  Group  and  select  time 
only,  we  should  be  able  to  see  the  result  like  Figure  28. 
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14.  Further  data  analysis  can  be  done  under  Report,  such  as  generate  a  movie  and  export 
data.  We  can  add  a  Domain  Probe  Point  under  Definition  to  analyze  the  temperature  rise 
as  a  function  of  time  at  a  fixed  point: 
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Click  on  Probe  ID  Plot  Group  under  Results,  we  can  plot  temperature  as  a  function  of 
time  at  the  selected  point. 
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